{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve the 1D grating Using Scattering Matrices\n",
    "\n",
    "Even in scattering matrix formalism, the eigensolver step is the same as the way we did it using the Monaharam formulation, so there's really nothing new here in terms of elucidating why there is no conservation of energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('D:\\\\RCWA\\\\')\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.linalg import cond\n",
    "import cmath;\n",
    "from numpy.linalg import solve as bslash\n",
    "from scipy.fftpack import fft, fftfreq, fftshift\n",
    "from TMM_functions import eigen_modes as em\n",
    "from TMM_functions import scatter_matrices as sm\n",
    "from RCWA_functions import redheffer_star as rs\n",
    "from RCWA_functions import rcwa_initial_conditions as ic\n",
    "from RCWA_functions import homogeneous_layer as hl\n",
    "from scipy import linalg as LA\n",
    "# Moharam et. al Formulation for stable and efficient implementation for RCWA\n",
    "plt.close(\"all\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4lFX2wPHvncmkB1KpARKK9CZFmgrYsCxiF3uXVVHXXte6rr+1rl3WXtYGVuwugiBNQHqHAAklpPcyydzfH3cmCSSZmSQzySScz/PwJJl55507L8mZO+eee6/SWiOEEKL1sLR0A4QQQjSMBG4hhGhlJHALIUQrI4FbCCFaGQncQgjRykjgFkKIVkYCtxBCtDISuIUQopWRwC2EEK1MkD9OGh8fr5OSkvxxaiGEaJNWrlyZqbVO8OZYvwTupKQkVqxY4Y9TCyFEm6SU2u3tsZIqEUKIVkYCtxBCtDISuIUQopXxS45bCHHksdvtpKWlUVpa2tJNCWihoaEkJiZis9kafQ4J3EIIn0hLSyMqKoqkpCSUUi3dnICktSYrK4u0tDSSk5MbfR5JlQghfKK0tJS4uDgJ2m4opYiLi2vypxIJ3EIIn5Gg7ZkvrpEEbnHE+nz2hyz4/feWboYQDSaBWxyRUrOLOXbdPVT87/GWborwoRdeeIH+/fszffp0TjzxRIYNG8Ynn3zC888/T3FxcUs3z2dkcFIckbKyMxmm8ulUuQ+ttXzEbyNeeeUVvv/+e9LT07n77rtZvXo1YGZzX3LJJYSHh7dwC31DArc4IpVn7ACgG+mUlFcQHtL40iwRGGbMmMHOnTs57bTT2Lp1K5GRkQwbNowrr7ySffv2MWnSJOLj4/n1119buqlNJoFbHJEqM1MAaKdK2J99kPDOXVu4RW1L0j3f+uW8u548vd77XnvtNX744QcWLlzI+vXrefrpp5k7dy4Azz33HL/++ivx8fF+aVdzkxy3OCJZ83ZVfV98YHvLNUSIRpAetzgiBRfsqfrenrkTOL7lGtMGuesZi6aTHrc4IkUUpbLFkQiAzklp4dYIf4uKiqKgoKClm+EzErjFESm6dC9bdDfSdTRBeV4vgyxaqeuuu45TTz2VSZMmtXRTfEJSJeLIU2kntiKdPXo0nXUWnWukTUTrtmvXLgAmTpzIxIkTq26fOXMmM2fObJlG+YH0uMWRJy8VKw52647s0R2JKklr6RYJ0SBeB26llFUp9adSaq4/GySE3+XsAiBVd2CPowNR9gywy1KkovVoSI/7FmCTvxoiRHPR2WYwcrejI7t1RyxoyJV0iWg9vArcSqlE4HTgDf82Rwj/s2emUKaDOEAMe3QHc6OzFy5Ea+Btj/t54C7AUd8BSqnrlFIrlFIrMjIyfNI4IfyhInMHaToBlIU9uqO5UUoCRSviMXArpc4ADmqtV7o7Tms9S2s9Ums9MiEhwWcNFMLXVM4uduuOdI0OI5N2FOsQyJbALVoPb3rc44GpSqldwMfAZKXUB35tlRD+ojW2/N3s0R3oFhNOkMXCbt0BhwRu4fTOO+9w0003eTxm3759VT9fc801bNy40d9Nq+IxcGut79VaJ2qtk4ALgXla60v83jIh/KE4i6CKIvbojsRGBBMdbmOP7iiBWzTI4YH7jTfeYMCAAc32/FLHLY4szkHI3boDMRE2YsKDTWVJ7i7QukWbJppu2rRpjBgxgoEDBzJr1iwAIiMjuf/++xk6dChjxowhPT0dgG+++YZjjjmG4cOHc+KJJ1bd7lJQUEBycjJ2ux2A/Px8kpKS+Oyzz1ixYgUXX3wxw4YNo6SkhIkTJ7JixQoAfvjhB44++miGDh3KCSec4JfX2aCZk1rr+cB8v7REiObg7Fnv0R0ZEB5MTHgwe7I6YKksg4ID0K5zCzewjfj+Hjiwzrfn7DQYTn3S7SFvvfUWsbGxlJSUMGrUKM455xyKiooYM2YM//jHP7jrrrv4z3/+wwMPPMCECRNYunQpSineeOMN/vWvf/HMM89UnSsqKoqJEyfy7bffMm3aND7++GPOOecczjvvPF5++WWefvppRo4cecjzZ2RkcO211/Lbb7+RnJxMdna2b6+Bk/S4xZHFWT2SqhOIDnelSjoccp9ovV544YWqnnVqairbtm0jODiYM844A4ARI0ZUTYtPS0vjlFNOYfDgwTz11FNs2LCh1vmuueYa3n77bQDefvttrrzySrfPv3TpUo477jiSk5MBiI2N9eGrqyZrlYgjS84u8qxxlBJCTLhJlSyrWcvdY1yLNq/N8NAz9of58+fzyy+/sGTJEsLDw5k4cSKlpaXYbLaqremsVisVFRWAWb/ktttuY+rUqcyfP5+HH3641jnHjx/Prl27WLBgAZWVlQwaNMhtG5prGzzpcYsjS3YKB6wmHRITEUxMRDB7dQIOLFIS2Mrl5eURExNDeHg4mzdvZunSpR6P79rV7Hz07rvv1nvcZZddxvTp0w/pbde3TOzYsWNZsGABKSnmd0lSJUL4Qs4uUpWZdBMTHkxMuA07QeQHd5TZk63clClTqKioYMiQITz44IOMGTPG7fEPP/ww5513Hscee6zbLc0uvvhicnJymD59etVtV1xxBTNmzKganHRJSEhg1qxZnH322QwdOpQLLrig6S+sDkr7YSR95MiR2jXCKkTAsJfAPzrxRtCFPF44lQV3TmTZzmzumrOWn2Ofpk+sFa75paVb2Wpt2rSJ/v37t3QzfG727Nl89dVXvP/++z47Z13XSim1Ums9sp6HHEJy3OLIkWM2TNhabnpXMc46boB9lk70yV7eYk0TgWnmzJl8//33fPfddy3dlENI4BZHDmcqZLs9niCLIiokiNiIYABTWVKcCWUFEBLVgo0UgeTFF19s6SbUSXLc4sjhLPfbrTsSHW4qDaLDTeDeUSmrBPqCP1KvbY0vrpEEbnHkyNlFpS2CLNpVBewYZ6pkc1lc1TGicUJDQ8nKypLg7YbWmqysLEJDQ5t0HkmViCNHdgqlkd2gQBHrDNztw2woBRtLYyEEKQlsgsTERNLS0pBlnd0LDQ0lMTGxSeeQwC2OHDkpFIT1AKgalAyyWmgXaiOvJAJHaAwWmT3ZaDabrWrGoPAvSZWII4PDATm7yQ7pApgabhdXuqS8XQ9JlYhWQQK3ODIU7IfKMg4GVc+adHHlu4vDEyVVIloFCdziyODsSe+lE1DdywaqSgJzwxIhLxUqK5q9eUI0hARucWSoKgU02+rVTJW48t1Zts7gqID8tOZvnxANIIFbHBmyU0BZSbGbsr/oGj1uVxB3LT4l6RIR6CRwiyNDzi5on0hmiQOoTo/U/H6P7lR9rBABTAK3ODLkpEBMEjlF5UD1gKT53vS+0yragTVYNlQQAU8Ctzgy5OyC2GRyis3+gTF1pEqySyohWkoCReCTwC3avtJ8KM6isn0P8kvtKGVmTLq4etw5RXaISZIctwh4ErhF2+dMfRRHdEdraBdqI8ha/avvynHnFJdDbLLpcct6GyKASeAWbZ8z9ZETataHqJkmMT+7ArcdYpKhLB9Kcpq1iUI0hARu0fY5Ux9ZNufkmxoVJVCdKsktLkfH9DjkMUIEIgncou3L2QVhsWRWhAGHTr4BCAmyEh5spcKhKYro7nyMBG4RuCRwi7avVimgrdYhVekSW+fqxwgRoCRwi7YvO8VZCmgC9+E9boCYCBPMs+1BENkJsnc1ZwuFaBAJ3KJtq7RDXhrEVNdwx0bUEbhdtdw1K0uECFASuEXblpcKuhJiksgt9pwqyS0uN7XckioRAUwCt2jbXD3n2GSyi9ykSg6ZhJMM+fvAXtpcrRSiQSRwi7bNVdYXk0Ru1XT32oE7+vAeNxpy9zRTI4VoGI+BWykVqpRarpRao5TaoJR6pDkaJoRP5KSANQSiulQPTkbUlSpxDk66ctyuxwoRgLzZLLgMmKy1LlRK2YBFSqnvtdZL/dw2IZouZxfE9ACLxUNVyWGzJ12PFSIAeQzcWmsNFDp/tDn/yUIOonXI3gUxSWitq1IlHgcnI+LBFiGzJ0XA8irHrZSyKqVWAweBn7XWy/zbLCF8QGtnjzuZgrIKKhyaiGArIUHWWodWlQMW2UEpKQkUAc2rwK21rtRaDwMSgdFKqUGHH6OUuk4ptUIptSIjI8PX7RSi4YqzoLzADEwWuXrbtdMkUJ33dpUMSkmgCGQNqirRWucC84Epddw3S2s9Ums9MiEhwUfNE6IJXKmO2GQz6EjdA5NQc4XAmoF7Fzgcfm6kEA3nTVVJglIq2vl9GHAisNnfDROiyVypjhj3090BwoOtBFstlNodlJRXmlRJRSkUpjdTY4Xwnjc97s7Ar0qptcAfmBz3XP82SwgfcKU6YnpUpUDqC9xKqeqdcKpquZF0iQhI3lSVrAWGN0NbhPCt7BSI6gy2MDPoSO1NFGqKjQjmYEEZOcXldHGVBGanQI9xzdFaIbwmMydF25W5FeJ6A9RYp6TuHre5zzVAaYfo7mCxmXMIEWAkcIu2SWvI2AIJ/YDqQce6VgZ0qS4JLAerDeL7QIYM54jAI4FbtE15aaYUsIMrcNc/+cblkPVKwAT9g5v8204hGkECt2ibXD3lDgMAqna/qW9wEiA2wjU4aa9+bO5uKC/yXzuFaAQJ3KJtcvWUEw7tcXudKoGq3joZW/zTRiEaSQK3aJsOboLIjhAeC+B2EwWX2qmS/tXnEiKASOAWbVPGpqreNuBxAg7UkSqJTTZLwmZI4BaBRQK3aHscDpPecOa3S8orKbU7CLZaCA+uvcCUS60et8UK8UfBQaksEYFFArdoe/L2gL24RkVJ9TolSql6H3bIhsEuHfpLSaAIOBK4Rdvj6iE7c9TepEnM/c4JOM5ZloAJ/nmpUJrv+3YK0UgSuEXbc3Cj+ZrQF3BuAIz7gUmAdqE2LAoKyiqwVzpXBXQNUEpliQggErhF25OxGaK6QFg04N2sSQCLRVXluauWd60qCZQBShE4JHCLtufgJpObdvJmnRKXQ9YrAYhOgqAwGaAUAUUCt2hbHJVmYagagdtV3uduZUCXWFeP2zUJx2KBhKOkxy0CigRu0bbk7DIbINSo4c72Yrq7S3WqpOYA5QDpcYuAIoFbtC2uWY51pEq8CdwxNTdTcEnoBwX7oCTXd+0UogkkcIu2xZXScFaUQI1UST37TdYUE3HY4CRUvwlIPbcIEBK4RdtycDO07w4hUVU35TRgcDKmavZkjVSJK+0ia5aIACGBW7QtGZurS/icqsoBG5AqqVohEKB9N7BFSI9bBAwJ3KLtqKwwFSUJhwbu3Kr9Jr0fnMytmSqxWMybgfS4RYCQwC3ajuydUFl+yMCkvdJBQVkFFgVRoR73xq4xOGk/9I6E/hK4RcCQwC3ajozaFSU189sWS/0LTLnE1jU4CabHXXQQirN901YhmkACt2g7Dm4GFMRXV5TkNmDyDdSo4y46LHDLpgoigEjgFm1HxiaI6QHB4VU3ebPXZE2uKe95JXYcDl19h6xZIgKIBG7RdhzcXN0zdmpIKSCAzWohKjQIh4b80hp57nZdIaSdzKAUAUECt2gbKsoha1sdpYCuTYK9S5VAHZsGAyhlqlUkVSICgARu0TZk7wBHRdV2ZS7ebqJQU72VJR36SapEBAQJ3KJtcPWED6/hLnZtouB94K6zlhtMGqY4CwozGt9OIXxAArdoGzI2g7KYzX1rcKU7GpIqqS4JrKPHDdLrFi1OArdoGw5ugphksIUecnNDNlFwcVWW1CoJdKVhZIBStDCPgVsp1U0p9atSapNSaoNS6pbmaJgQDXLYrjcu1ZsoNCTHXc8knMiOEBpdvaelEC3Emx53BXC71ro/MAa4USk1wMNjhGg+FWVmunudgds1ONmAqpL6UiVKmeeQxaZEC/MYuLXW+7XWq5zfFwCbgK7+bpgQXsvcBrqy1sAk1JiA42Gj4Jpi6kuVQHVJoNa17xOimTQox62USgKGA8v80RghGsXVAz6sx+1waPJKnFUlYQ2v466VKnE9R2kuFKY3rq1C+IDXgVspFQnMAW7VWufXcf91SqkVSqkVGRlSLiWa0cFNoKwQ1/uQm/NL7Ti0WRUwyOp9H6XWTu81dZA1S0TL8+q3WSllwwTtD7XWn9d1jNZ6ltZ6pNZ6ZEJCgi/bKIR7BzeZoB0UcsjN1aWA3qdJah5fZ49bFpsSAcCbqhIFvAls0lo/6/8mCdFAGZtqTXWH6sHFhpQCwqGpEn14LjsyAcLjpJZbtChvetzjgUuByUqp1c5/p/m5XUJ4x14C2Sm1FpeCmru7e5/fBgi1WQm1WbBXaorKK2sfkNBfarlFi/K4JYjWehHgeQV6IVpC5lZAu+1xN6SG2yUmPJj9eaXkFJUTGXLYn0mHfrD2U1NZouRPQzQ/mTkpWreqNUrqqOFu4FrcNbmtLEnoB2X5kL+3wecVwhckcIvW7eAmsNggrletuxoz+cYlJqKeFQJBpr6LFieBW7RuGZshvg9YawfnqsHJBlaVgJsVAqG6JFAGKEULkcAtWreDm+qcMQnVqZLYRqRKYuvaTMElPBYiOkiPW7QYCdyi9Sovgtzdda5RAk1MldS3mYKLbKogWpAEbtF6uaa619PjbswmCi5uUyVQXRLocDT43EI0lQRu0Xq5UhUd6l6ssqrH3YBNFFxcj6kzVQKml28vgrzUBp9biKaSwC1ar8wtYA2GmKRadzkcutFT3gHiIsz0+XoDt6uXn7m1wecWoqkkcIvWK3M7xPYEa+15ZHkldiocmnahQYQEWRt86vhIE7gzC8vqOaCPsw3bGnxuIZpKArdovbK211oR0CXDGXDjo0LqvN+T+CjTS88oqCdwh8eZ3XCytjfq/EI0hQRu0TpVVphdb+oJ3JnOgOvqOTdUbHgwSpmqEntlHQOQSpnnzpIet2h+ErhF65S7Gxz26pTFYVw97oRGBu4gq8V9LTeY586UHrdofhK4ReuUtcN8ra/HXWiCbXxkwwcmXVy99XrTJXG9oGAflBU2+jmEaAwJ3KJ1cqUo4urucbsGFRMameOu+dh6Byhdz529o9HPIURjSOAWrVPmNgiLgYi4Ou/OaGKO2zzWwwClVJaIFiKBW7RObipKoLqX3LTA7epx15Pjju0JqOq0jRDNRAK3aJ2yttebJoEagbsJqZJ4T6kSWxi07yaVJaLZSeAWrU9ZARTsh3g3Pe4C3w1O1hu4wbRBUiWimUngFq2Ph4oSh0OTVdT0VInHwUlXG7J2mG3MhGgmErhF6+OarVhPqiSvxI69UhMVGkSoreHT3V1cvXVX771OcX2gvAAK0xv9PEI0lARu0fpkbgOUc3CwjrubOPnGxfX4DE+pkqo2CdE8JHCL1idrO0R3A1tonXdn+KCiBMyqgmbaezkVdU17h+p0jaxZIpqRBG7R+mRt81BR4hyYjGr8wCSYae8x4cFo7Wbae7tECAqTwC2alQRu0bpobQYD61mjBJq+wFRNVZNw6kuXWCxm6rukSkQzksAtWpeCA1Be6NXkm6bmuKFmZYm7Acpe0uMWzUoCt2hdqtYoqT9wV013b8LkG5eqWu76pr2DSdvk7IIKN8FdCB+SwC1aF1dKwl2qxEeDkzXP4b6ypA/oShO8hWgGErhF65K1wwwGRnWp9xBfLOnq4l2PWypLRPOSwC1al6xtJlBa6v/V9W2PO/iQc9apKnDLAKVoHhK4ReuSuc3tGiVaa5+sxe3i1eBkWDREJEhliWg2HgO3UuotpdRBpdT65miQEPWqKDdblrkZmKya7h7StOnuLl4tNAXVa5YI0Qy86XG/A0zxczuE8CwnBbTD78u51uTqcde7mYKLbBwsmpHHwK21/g3Iboa2COFeVUWJu1JA3w1Mgpn2DpDtbto7mMqSogwoyfXJ8wrhjuS4RetRtSqgf3e+qclmtRATbjPT3ovdTcJxDVBKukT4n88Ct1LqOqXUCqXUioyMDF+dVohqWdsgogOEtq/3EF8H7prn8ri8q6uNQviZzwK31nqW1nqk1npkQkKCr04rRLXM7W4n3kB1LtoXFSUuXm2oEJMEyiqVJaJZSKpEtB5Z2826IG74tcftLnAHBUNMD5mEI5qFN+WAHwFLgL5KqTSl1NX+b5YQhynJgeJMtxUl4NtZky5V0949Vpb0kcAtmkWQpwO01tOboyFCuJXpDIgeUiW+Lgc05/Ji9iSYtqX8Bg6H25mdQjSV/HaJ1sGLihKoXlPEF0u6ulSnSjys/hfXCypKIH+vz55biLpI4BatQ9Y2M/gXk1TvIWa6uytV0syDkyCVJaLZSOAWrUPmNhO0rbZ6D8kvqaC80kFkSBBhwU2f7u6S4HWO27VxsOS5hX9J4Batg4ftyqDmJsG+G5g05/MyVRLVCYIjZYBS+J0EbhH4HA7I3uE5v+2HUkCAOOcbQXZRGZUOXf+BSsmaJaJZSOAWgS8vFSpKWyxw26wWosNtONzt9u4S11tSJcLvJHCLwJflZSlg1V6Tvk2VQAOWd43vY95o7CU+b4MQLhK4ReCrKgX0LsedEBnq8yYkNGRdbjRk7/R5G4RwkcAtAl/mNgiOgsgO7g9zLenqjx631yWBrsoSyXML/5HALQJf1nazBrdSbg/zV47bnNO8GXhdEiiVJcKPJHCLwJe13WOaBPwduL0sCQyJNDvQS+AWfiSBWwS28mIz2OehogSqg6ovp7u7VOW4PfW4wUx9l1SJ8CMJ3CKwuQb53GxXBma6e4Yfq0qq9p70lOMGU1mStQ20m5pvIZpAArcIbK7JLB5SJfmlZrp7RLCV8GCPi142mNepEjBtLc2D4iyft0MIkMAtAp1rMou3Gyj4cDnXmrxe2hWkskT4nQRuEdiytkO7rhAc4fawqsk3fshvA8RFmPNmFXqY9g7VaR0ZoBR+IoFbgKOyRZ5Wa81PGw6Q5a4Xm7WtQQOTvl5gyiU4yEL7MDPtPcfdbu8A0T3AYnO7Zsm29AJW7Mr2cSsboIX+z4VvSOA+kh3cDJ9eBk90NTu3NLMPlu7muvdX8sg3G+s+QGuTKvEqcPu3x23O7WW6xGKF2J71rlnicGgufXM5572+hM0H8n3dTM8WPmv+z396AIokD98aSeA+EmXtgDnXwitjYPs8CI+FOddA4cFma0J5hYNX5+8A4LdtGTjqSj8UZUJZnsc1SsA/u7sfrmpDhQIvBijj699/cvOBAg7kl6I1vPzrDl820bOU32DeY9A+EZa8DP8eAvMeh5Lc5m2HaBIJ3EeS3D3w1U3w0ijY9A2MvwVuXQsXz4bSfJhzdbN9hP7izzT25ZWaZhXb2bi/jp6nlxUl0Fw9bi+nvYP5lJC9Eyrtte5avCOz6vu5a/exI6PQZ210qyAdZl9t2nbdfLhhKfQ5CX57ygTw356CsoLmaYtoEgncR4L8/fDtHfDC0bD2Uxh9HdyyBk56xPS2Ow6A0582vbHfnvJ7cyoqHVU9zY7tTDBcsqOOj+x7V5mvnQZ5PGfABe5Og8Fhh/QNte5yvdaO7ULQGl5pjl63oxI+v8YE5vPeNTM8E/rCee/AjEXQY4Lpef97KCx+0Ux8EgFLAndbl7nd9LBXvg1HXwo3/wmnPglRHQ89btjFMHQ6zH8Sds73a5O+WbuPPdnFJMdHcMfJfYFDe6FV0pZDdHezs4wHGa5Zk36YfONSNQnHm9mT3Uabr2l/HHJzRaWDZSlmUPKFC4djtSi+XL2X1Gw/B8oF/zJvzKc9Zd6oa+o0GKb/F66ZB52Hmdz3GydChRevU7QICdxtmdYw91awWODG5XDGc9C+a93HKgWnPwPxR5l8d8EBvzTJ4dC8NM/kfv86sRfH9kkAYHlKNvZKx6FtT10O3Y7x6rz+Lgc053YuNOVNj7t9N4jqDKnLDrl53d48CssqSI6P4JiecZw5rAuVDs2rC/zY697xKyz4P/PGPPyS+o9LHAGXfm565Ac3wKLn/dcm0SQSuNuyNR/DroVw4iMeJ7AAplb6/HehvMgE78oKnzfphw0H2JFRRNfoMM4a3pVO7UPpmRBBUXkla9Pyqg/MS4WC/ZA42uM5tdY19ptshsFJb2ZPKgWJo2oF7sXONMnYXnEA3DCxN0rB7BVp7M/zw+YLBQfg82tNWuT0ZzyusAjAwGkw6FxY+LTs5hOgJHB7Y88y+PhimHsbZGxt6dZ4pygLfrzP9FiPvtz7x3Xob/7Ady00vTQf0lrzorO3PWNiL2xW8+s3zhnEltRMl6QuN1+7eQ7cBWUVlFc4CA+2EhHi++nuLvENWWgKzLXP3XPIpxdXftv1mnt3iOS0wZ0pr3Qw6zcfb75QWWEGI8uLTC/awySmQ5zyBNjCzCe21rDmitawcwF8ejl8fbO57m2YBG530jfCR9PhrZNhzxL48wN4eRR8cC7smBfYv9A/Pwhl+XDG8yZV0hDDLoJhl5iByu3/81mT5m0+yKb9+XSICuG8EYlVt4/rFQ9U90YBE7ht4dDRi4HJZkiT1Dy/V4OTUJ3mcb4JlVVU8odz0s2YnnFVh900ydSpf7R8j3f5c28teBJ2LzJvxB36NeyxUR3NJ7VdC2HNR75rk6/ZS83f5WsT4L2pJo+/5mN4cQT8cK8pKW2DJHDXJXcPfPFXeHUc7FoEkx+AW9fB3zbAxPtg/xp4/yxTB73yHb/vL2ivdFBRM//rScpCWP0hjLu59kCUt057ChL6wefXmaqUJqrZ277uuJ6E2qxV97mC2IrdOZTaneWIacuh6wiweu5B+3vWpItrt/esovK6684P13kIWIPNawH+3JNLWYWDfp2iDnmT6d+5HScN6Eip3cGbi1J809jt/4PfnjZvwMMuatw5jr4cuo2BH+9v0ESdUnsl2t+dmoJ0+PUJeG4gfHWj6URNfQlu2wQ3r4Ih58Oy1+Dfw2D+/7W5Mkf/fa5sjQozYOEzsOJNQMG4m2DCbaZkDsxHzYl3w4RbYf3nsPRl+OYW+OURGHmlOTYk0qdNyiws4/zXl1BYWsHLFx/NqKRY9w+oKDMfb2OS4Lg7PZ6/uLyCnzemU1hWQZndQVmFg1J7JWUVDqISHuTarKtIf+MiOt38C7Ygq8fz1ef37VmsTs0lNiKYi47pfsh9sRHB9O9GP41gAAAgAElEQVTcjk3781m1O4dx3cNg/1pznb3QHKWAACFBVtqFBpFfWkFOcTlxnp4vKAS6DK/qcS/ebnp/rvx2TTdN6s3PG9N5f8kuZhzfk+jwxr8JHdyfSsTHV1EYmsRb6mqsP2wmJMhKqM1CSJCFEJuVDlEhTO7XAeUu522xwF+eN73Znx+Eaa+4fV6tNR8s3c1j327ihH4deOmio7FavMipN0R2iqmQWT/b1MgfNQXG/BWSj6vO37dPhDNfNh2XeY/B/Cdg+Szz9zDySvP/0spJ4HZZ9jr871GwF5uR9+Pvqb8CIygEhk2HoRfC7t/Jn/8CUQufpSBlBe2u/NyrXqI3Su2VXPfeCnZmFAFw0X+W8uiZg5g+unv9D1r0nJmxd8nnEBzu8TnunrOOb9bsq/f+dOsFPOp4l59++YaTp0xr8GtweWGemUxz9YTkOpddHdcrjk3781m8I4txQRmgK72uKKleh9v/f5DxUSHkl1aQWehF4AaTo1/2OlSUVaWCXKmhmoZ2i+a4oxL4bWsGb/2+i9tOOqrRbVw25zn+UpHL2UX3sGVJ/Z+WHp82iEvG9HB/sg79zUSthc+YqpTkY+s8rLzCwUNfb+Cj5Sa3/P36Azzx3SYePKORn/jqUnCA4v9MIagsj8ohlxJ27I3uB90T+sIFH0DaSvjlIfjhbjNbdNrLJtC3YpIq0drULn9/F3QfCzcsg6kv1h+0a9ifX8rty6MYuuVy7rVfTbu0BWR8dotPct9aa+6avZZVe3Lp0j6US8Z0x16puffzdTz01fpDS+dcMraaP7DB50HvEzw+x4KtGXyzZh+hNgvTR3fjinFJXH98T24+oQ93ntKXB88YQPjoyyjUoVT+8U7D0jU1LNuZxfKUbNqFBnHZ2LoDhWuwbvGOzOqBycRRXp0/s2p3d/8Hbq93e3dJHA2V5ZTsWcXq1FwsCkYn1/2paeZkk+t+5/cU8ktrz7j0xu7MAoZlfM0SxwDOP30KD5zenztP6cvMyb257rieXD62B6cP6QzA//2wmYMFpZ5Petyd5hPc3FtNTvkwmYVlXPLGMj5avofgIAs3TOyFzap4c1EKHy7b3ajXUUt5EZn/OQtdnMNZJQ8yctXJvLDaQXG5F5VPiSPg8m/g0i/AFgofngfbfvFNu1pIm+lx78go5OvV+7AoRUSIlbBgK+HORfVdX7vFhtEhKrT6QVqbj1ILnzETUKa+aBYI8qCg1M7rC3byxqKdlNod2KyK9Z2m8Vp6OjM2f8C+H5LpcuodTXo9//7fNr5es4+IYCtvXjGK/p3bMSQxmge+WM+7S3azNb2QVy4+mpiI4OrXMvdvphLglCc8nr/UXsmDX64H4NYTj2LG8XX3XCodScxddzwn23/l51VbOXVUAwe5gJd+NbntK8cnExVqq/OY0cmxWC2KNWl5VLRbRlBcn+oUlQf+Xou7Jq93e3dxVsXsXbeACsdghia2p31Y3ddgVFIsxyTHsiwlm/eX7ObGSZ4X1zrcvO8+5UqVwYJuN3D1hOQ6j9FaU1RWwfwtGTw+dxMvTB/u/qS2MDMH4P2zzCe6SfdW3bVhXx7XvbeSvbkldGwXwqxLRzK0WzRJ8RHcNXstf/9qAz1iI5jQp/anDK85Kkl98xK65G3iWvvtBHcbRtGeXJ79eSsfLtvN7Sf35ZyjE92nZZSCXpPhyh/g/TPh4+lw/nvQ99RDDtuZUUhWUTlFZRWUlFdSXF5JcXkFxeWVFJVXEh5sZfro7vX+HzYXrwK3UmoK8G/ACryhtX7Sr61qoM9XpXH/F+spsbtfZ8NqUfz9jAFcNrYHCkzebvGLZhDGi+oLe6WDj/9I5fmft5JVZAbEThvcibtO6UdiTBg3/zeU77amM2XZ4+xp34Pu485r1Ov5avVenv9lGxYFL0wfTv/O7QA4f2Q3eiVEcv37K1myM4upLy/iP5eNpF+ndrD6v6aC4C//hsgOHp/jxXnb2JNdTL9OUfX+gYO5ZsGjryRs8Y+kzHsbPfJJ93nRw6xOzWXhtkwigq1cOT6p3uOiQm0MSWzPn3tycKQug/6ne/0cGQWuvSb9OzhpnqMBsyfBzPqM7kHFrqXAYMb1dh/AZk7uw7I3l/HmohSuHJ/UoN18DhaU0mn7J2SrSMacflm9xymleOzMQZz03AK+XrOPc0ckctxRCe5P3msyDD4fFj0Lg86BhKP4du1+7vhsDSX2SoZ1i2bWpSPo0M50jM4f2Y2dGUW8tmAHf/1wJV/cMJ7eHRo3/rPlg7/RN30eD1dcxvjTLuHqCcks2ZHFE99tYt3ePO6avZa3FqVw32n9Pb+OiDjT+37/bPjkEjj3bRgwlVJ7JQ98uZ7ZK9M8tufDZbt5+aKjGZIY3ajX4wseUyVKKSvwMnAqMACYrpTyYeKqWkNHokvKK7l79lpu+9T88kwZ2ImbJvXmqvHJXDiqG1OHduHE/h0Z1yuOIYntqXRoHvp6A/fOWUvFd3eZoD3qWo9BW2vNLxvTmfL8bzz45Xqyiso5uns0c/46llcuHkFSfARBVgv/vmgE3/Z6iLWOniT8dCMp6xY1+Bqs2pPDnbPXAnD/6QM4of+hU9NH9Ijhm5njGdy1PanZJZz9ymLmrdpopil3GwPD6/+DddmWXlBVM/yPswZV1VPXZ9Kkk9lMMpOKvuO3rRkNej2uWZKXjk3yOOA2rlccSeoAwWU50M27NAk03+CkeY4GzJ506TaahNy1gK5KCdVnfO84hnePJruonP8ua1gt8ie/ruREtYKV0afQu7P7N4huseHcfIJZvOvBr9ZXV/O4c8oTYAtHz72VZ37czI3/XUWJvZJzjk7k4+vGVAVtl7tO6cspAztSUFrB1e/+QXaRFxOXDrP6i2fpu/Nd3qk4mc4n31rVyRjbK46vbhzP8xcMo2t0GJsPFHDZW8u59M1lbKprwbKawmLgsi+hy9Hw2RXk/vExF85ayuyVaYTaLIzoEcOxfeI5ZWBHzhrelYuO6c61xyZz8wl9GNilHanZJZz76hLeXbzL/9Uz9fDm7Xw0sF1rvRNAKfUxcCZQzyLKjVNcZufbp6+irNcUxp1wJj0T3L87bz9YyI0frmJLegEhQRYemTqQC0Z1c9sb/Gr1Xu6Zs5qBqx8lKOgXikdcT/hp/+d2NtmWAwU8Nncji5wVAUlx4dw9pR9TBnWq9Vw2q4VnLxnL3e/8kztSbyByziWkhP9Icq++Xl2D1OxirntvBeUVDi46pjtX1dND7dw+jM9mjOXuOWv5avU+cr64k8qgfCxnPIfy8KnB4dDc/8V67JWa6aO7M6KH53REqM1Kep/pHL/tCT74+TuO73uFV69nyY4sftmUTqjNwjXH1t+rdxnXK54DC5wrAno5MAnNHbgbsLSrU0nHkcSt+4we1ixGerjeSilmTu7NVe+s4JX5O5g2vKtXryu/1I595YfYVCVdT/irV+269tiefPnnXramF/Lyr9u5/WQPv6eRCZRPfpjg727lwPY3saiJ3Hdaf66ekFzn353FonjugmGc//oS1u/NZ8b7K3n/mtGEeFmdtOznTxmx+jHmOYZRNOlRbjwsnWexKKYN78qUQZ14Z/EuXp63nYXbMjn9hYVMH92d20/uS2xEPZ2F0PZw6ecUvnU2Ud/+laTyGWREn8R/LhvJgC7t6m3TDRN78cR3m3hvyW4e+noDy1KyePKcIbSrJwXoL94MTnYFUmv8nOa87RBKqeuUUiuUUisyMhrWKwNYuG4bI8qWM33Tjcx+/jYufO13vvxzb509ga9W72XqS4vYkl5Az/gIvrxxPBeO7u7xI/yZQzrx+4CvuTToF16t+AsnrDuZdXvrfnfOKSrn71+t57QXFrJoeybtw2z8/YwB/PS34zl1cOd6nyskyMqTl5/Iy53+QYguxf7+eaTs9bzuR0GpnWveXUFmYTkTesfzyNSBbl9PqM3K8xcM483hOzjHuoiX7H/hvt8rPA4gzl6ZxvJd2cRHBnPPFO/z1cPPuJZiHcKgA1+yak+OV6/njs/WADDj+F5eBZ8RPWIYad1Ovg4nJ9xzoIdDd3f351rcLgkNzXEDf2J6tmcl7CMs2HPQmtS3A+N6xZFdVM79X6zzqlf34ZLdTNP/Y3PwQAYM8e7Tis1q4YmzBgPw2oIdbD/ovtY5q7CMC5b3ZqmjP4/Y3uWTs2O55tiebn9Pw4ODeOOyUXRsF8LyXdnc9/l6r17Pot8XMGDRzWzV3dgw7nluPKF/vceG2qzMOL4XC+6axBXjklBK8eGyPUx86lfeWpRS90A+MGd9HuP33sCyyn48G/wqPx6/x23Qdj3Xo2cO4qWLhhMZEsR36w7wlxcXsX5vntvH+Zo3gbuu/5VaV15rPUtrPVJrPTIhwUOeqQ4nj+hP3mW/sLb9ZO6yfcL1e+/l4U8WcswT/+Phrzew+UA+pfZK7v18Lbd8vJri8kqmDu3C1zMnVOWA3XJUwpc3ELvlY4rG3M7/uvyV/fllnPvaYr5avbfqMHulg3d+T2Hi0/N5b4kZEb9sbA/m3zGRqyYkExzk+ZKF2qw8dM15vBz/ID11KvvfmM7ujPr/YysqHcz86E+2pBfQKyGCly8+2mP6AkBlbOaE7f8kO34Ur3MuHy1P5fr3V9Y70p5VWMYT328C4IHTB9A+3PteQrv2sWzvcDJTrYt5e946j8c/+s1G9uaWMCSxvdeDbKE2K+NCtvOnozdLUzy/OQCm/rzCQZjNv9PdXRo8exL4OTOOYh3C8WHeTa5RSvHUeUOJCgnixw3pfL5qr9vjS+2VrF44l56WA+gRV3jdLoCRSbFcOKob9krzSay+oLonq5hzX1vCn2n5/DPsdoLDIhm1/FYznd6DTu1DefPyUYTZrMxZlcYr890vqPXbyvUk/3QVRYQyf8SL3DRlmFevJTYimIenDuSHW47l2D7x5JdW8OjcjUx5/jfmb6neJKSi0sHjczdy+2dryKsI5pfhL0LPSUT+eAuseNur5zpjSBfmzpzAgM7t2J1VzNmvLOb9Jc2XOlGenkgpNRZ4WGt9ivPnewG01v+s7zEjR47UK1asaFyLtKZ0ySxsP99Ptorm+pIbWaVNTWt0uI3cYjvBQRYe/stApo92nxqpkrHFrDOyexFMegCOv9NZd7qej5abDxPXH9+Tcb3ieXzuRrYdNAvbj+8dx9/PGEjfTlGNeinF5RW8/9IjXJ//Ap9ZpvBDt9spd2jKKhyUVzicXyspKqvkQH4pMeE2vrxxPD3ivFhToqwAZk2C0jyYsZCV2cFc/e4KcovtDOsWzZuXj6xVZ3z7p2uYsyqNCb3jef/q0Q0aZATI3rKI2I9O51771Vx9yyP07lD3dfl5YzrXvreC4CAL3908od7jainNQz/Zg+fs55Az6m88Ns3zdPeUzCImPT2fbrFhLLxrckNeTqPszS1h/JPz6NguhGX3nejVY056dgGP5NzD0A5WImZ6P+7x2YpU7py91gTwvx1Hl+iwOo/777I9RMy9nslBa4i8dzvKi/r9mnKLyznhmQVkFZXz1LlDOG9kt0PuX5eWx5XvLCezsJwBndvxzpWj6JC5FN6bZkpPz57l1eJVP244wIwPVqI19IgLJ9hqITjI/AsJshAcZCVSlTEj5RZ6qzQ+Gvg6V503rcG/p2A+ic3bfJDH5m5kV5ZZMndyvw7MnNybZ3/eysJtmQRZFI+eOchMCLOXmm38tv0II66AEx7yqqqp1F7J499u5IOlZjzi9CGdefLswfVWT7mjlFqptR7pzbHe9Lj/APoopZKVUsHAhcDXDW6Vt5QidNz1WK/9mYT2EcwJfZxXei4hKsRKbrGdpLhwvrhhHBcd4zk1QnmxmVTz6nhIX29mUx1vZhMGB5mPiY+eORCrRfH6gp1c/tZyth0spHtsOLMuHcEHVx/T6KAN5mPixTc+xFfh53Ce4weG73iJhdsyWJ6SzerUXDbtz2dHRhEH8kuJCgni9UtHehe0tTYzNrN3wLlvQVQnRvSIZfaMcXSNDmN1ai7nvLqY3VnVvaElO7KYsyqN4CALj00b1Kg/htijxnMgtCcXWn/l9QV1L4iUVVjGvZ+bwdW7p/TzPmgD7F2JQrNK96l7fe46NGd+GyDOmTPNKvRu2vvBglK2HSxkrTqK8OyNXvVQXc4dkchJAzpSUFbBnbPX1Pl8lQ7NRwv+ZIplOZnJ0xoctAGiw4O5/3STinjiu02HDCLO33KQC2YtIbOwnGP7xPPJ9c5ByJ4TYdL9sO5Ts9a7F04Z2IkHTh+A1aLYnVXMtoOFbNiXz597clm6M5vlW9O4aOfdDFI7+arXo40O2mA+tZzQvyM//e147jutH1EhQczbfJCzXlnMwm2ZxEUE899rx1TP4rWFmsk6Y2+CVe+bNezXfOxxTkaozcrj0wbzwvThRARb+Xbtfv7c4/9t4Dx+ttRaVyilbgJ+xJQDvqW1rr2th691GQ7XLUB9dSOnbX6RU/ruZOmQxxjWp4d3H4m3/gTf3QG5u2HoRXDSoxB5aApHKcVlY5Po0yGKG/+7ijJ7JTNP6MOV45O8HkDxJDIkiFNufZ29nwRx045PmDqkI3uPvptgm9VMP3b2OOIjQ7z/qP/HG7B+jukV1JjJ1rtDJF/cMI4r3v6DjfvzOefVxbx1xSj6dori/i9NeuPGib1Jjm/AKnE1KVMaOPS3B/n76kXsO+moQ3qBWpuP25mF5YzpGcuV45Iadv7U5WgUW4KOIiOjiPT8UjoeVqlwuOZaYMol1GYlKjSIgtIKckvs9Q9+OblWAyzsMAKV+SXs+xOSJnj1XEop/nn2YFbuzuH37Vm8t2QXV4w/NPf//fr9jMr7mRBbBd1O9G5Qsi5nDe/K7JVpLHaW2T193lA+W5HKPZ+vo9KhOXt4V548Z8ihqcJjbzfL1n5/t/l77eKhHhwzc/as4V3JK7FTXvXJs5KK0kL6zruW6IMb2TnhKS484ZpGB+2agoMsXHdcL84ansgzP23hkxWpDOjcjlmXjaTr4Z9ggoLhlH+YGdFz/wZfXG8WsDr9WUhwP5N16tAuDOrSjkXbMz2XJPqAV5FCa/0d8J2f21JbWLR5F1z6KtafH2T8wXWQfj50HGh27YjtWXvCTN5e+OEe2PS12RTg8rn1TtN1GdsrjoV3TUJjAq2vhQbb6Hrxa/BdJN1XzKJ7+2A4+XHv1kY+XNpKs+rZUVNgfO21PDq0C+WT68dww4erWLgtkwtnLWVi3wR2ZhTRMyGCGRN7Num1xI69lPKFj3GumsebiyYdMqX5y9V7+WHDASJDgnjq3KFYGrpORepyVMeBDAxNZP6WDJbsyGLacPczWJtjHe7DJUSGUFBaQWZhmdeBO67feFiEmRXqZeAG87qeOGswMz5YyZM/bObYoxLo5ay40lrz6q/bed46j8z2g4nvMrjRr0kpxePTBjHl+YXMXpmGQ+uq3PoNE3tx5yl9awdSi8WkSV471qQZrv/NlNp5EBsRfOh1KyuE/86AjOVw9ix6DTm/0a+jPglRITx5zhDuOKUvMeHB7ifrdBoMV/0Ef74HPz9kFpsbfwscd4eZjFSTvRQyNkH6BnoeWE9PexHwos/bf7jAnzmpFIy9wUx//u4OM3NLOytNgsLMWgqdBpnlP+3FZkU0RwVMftAsMhPk3aQMvw9sWSxmeU2LFZa8ZAZLp/yzYcG7OBs+uxzadYZpr9Zbex4VauPNy0dxz5y1fP7nXr5bZ6pa/jFtcNM/SYTFUNz7DM7c+j2Tlm/lpkm9iYkIZl9uCX//ynwQ+/sZA+gW28CP7A6H2eZr0DmMax/H/C0ZLN6R6TFwZzZjRYlLfFQIOzOLyCwo46iO7lNBrvVJRvTrBZv6VE/nb4Apgzpx9vCufP7nXm77dA1zZowlyGph4bZMQg+soE/IXuwT7m7Ua6mpZ0IkN0zqxfO/bOPzVXtRCh6dOpBLxybV/6DwWLP5xltT4IsZcOFHDVtGuKzATEFPXQ5n/wcGn9vk1+GO12/wFovJdfc93UzUW/g0rPvMBO+iTJN6PbDerAvkike2cLP1m9aN65Q1QOAHbpduo+D6BWb1u4zN5qKlb4D0dbBpLqx6zxzX+ySzJGmsd+VkzUopOPVfoKyw7FXzH37qv7z7T3Y4zBKrhelw1Y8eB06Cgyw8c/5QOrUP5ZX5O7jomO51rkrXGNETroVtnzOpYjHvLRnIzMm9uWv2WgpKKzixfwfOG5no+SSHy9hs1g/vdgzjEupYn7u+hxU236xJl6rZkx4qS1Kzi9mTXUxUaBADu7Q3telbv2/UH/ZDUweyZGcWa1JzeW3BDm6a3IdX5+9getCvlFsjCB7qm4D314m9+HFDOimZhTx/wXCmDPK83yeJI83knO/vhMX/hgl/8+7JSvPhg3Ng70o4900YeFbTGu8PkQlw1mtmSYxvb4OvZ5rb23c3HcYBU02nsdNgs56LF0tm+ELrCdwuQSHQeaj556K12WWkOMukUfz8btckSpmeds2e92lPe+6lLHoGtv9s8m1dj/byqRR3TenHVROSqwbVfKL7WIrb9WJ67jyuWXwiITYLi7ZnEhsRzD/PHtK43GRa9Y43/WPa0T7MRlpOCanZxW577809OGmeK9j53O4n4SzZad54xvSMMx/Nu42G1R9A9k7vtpKroX2YjafOHcolby7j+V+2ER0ezIade3grZCkMnt6w3W3cCAmy8sUN47BXOhpWGTH6WrPZyP8eNZ+OPaWDSvPMtPP9q81O8wOmNqndfpd8LMz4HQ5uNAE6rOWmu0NbWR1QKZM+6DQosIO2i1Imxz3+FrP297d/Mz3qmirKzIYOqcvhjzfNovGDz4eRVzX46eIjQ3wy0FNFKcLGXMkIyzbiS1J48vvNAPxj2qDGpyxSl0N4HMT2xGpRjOlpPlF4qi5pzgWmXLyt5T58m7KqbdgO24fSWxP6xHP52B5UODQPfLmeM62/E6bKCR7d8N8Jd8wAbAPL2ZSCqS9AbC+YfRVs/tYMxBYcMJ2TmkpyTSnh/jVmoadAD9ouQcHQZViLB21ojT3utkIpszWUJcisTpi1w3yaKDhgNsktPixN0GGgWaEtQN6Y1NCLcPzyCNOt83i04jLOGt6VUwd3bvwJU5eZVILz9Y3rFc+PG9L5fXsWF4yqf/3xjGauKoHqNwl3C01prfnduUxC1frb8X0hpL15rY3cleaeU/uzcFsmOzMLuThoHvYOQ7B18W6Cit+FRMEF78MbJ8HHNV6fskBkR7PgVlRn84kje6c59rDV+YR3JHC3JKXMIKotDFa+Z1Yui+5uemZRnat/0aM6mQqZQNq5IyIO1X8q52/8kc9iruHhvwxs/LmKsswgz/BLqm6qXp87C611nZ8YtNbVa3E3Y4/bmzW5d2QUcbCgjLiIYI7q6Fx3x2IxYzWpfzT6ucOCrTx7wTCeevND+rEHRnuZT24uHfrD39ZB1k4odHZCCmp8zd1jdq654EM46uSWbm2rJYG7pSllFqr3YpuxQKNGXE7khjl8d2I2qgHT52tJcwayxOod3Xt3iCQhKoSMgjK2HyykTx3VG0XllZTaHYTaLER4sQaIr3izJrdrx/qxveIOfdNJHA3z/2lyvKHtG/X8w7pF88HwzbA+HAb5twqjUcJizOYFwm/aRo5btIykYyG2J+qP/zRt15+05SZlVGMCh1Kqqtf99E9b6lx/pebkG5/m8D2oGpysZ4XA1Oxi3l68C4Dxh6+/3W00oE0lRWMVZ6PWzzHrYod6sU6PaHMkcIvGs1hg3EzTY978bePPk7ocOg2ptUfmleOTiXQutHTOq0tIzS4+5P6WqCip+XxZRWW1pqH/vj2Tv7y0yEx4io/gtEGH5f27jjA530bUc1dZ8C+oKIGxNzb+HKJVk8Atmmb4ZWbQ7ZeHTO6yoSorTO+z2+hadw3rFs0XN4wjKS6cTfvzOfPl36sqNaDlAneozUpUSBD2Sk1eiXnNWmveXJTCZW8tJ7fYzqS+CXxx4/jaKzCGtoMOAxpdWUL2TrPkwfBLTT5ZHJEkcIumsQbBSY+YwcWV7zT88enrzYzXOgI3QJ+OUXx14wSOPyqB7KJyLnlzGe/8nnLYOtzNN/nGpWaeu9ReyR2freWxuRupdGhumNiLNy4fVf++hN1GQ9qK2iWg3vjlEbAGw6T7mtB60dpJ4BZNd9QUk++e/6SZDdcQrpSBmx1v2ofbeOuKUVx/fE8qHZqHv9nI3XPWsjfX7DjeHLu7H871nOv35XHBrKXMWZVGmM3KSxcN564p/dyvhdHtGDNLNGNzw540dTls/BLG32wqjcQRSwK3aDql4OTHoDgTfn++YY9NXQZRXaC9+2nyVovi3lP78+8LhxFqs/DpijRm/WYW5G/OyTcu8c5e/u2frmFNai5do8OY89dxnDGki+cHJzp3qGlIukRrs69oZEez9Kg4okngFr7RZbhZVH/Jy2aFRm+lLa83TVKXM4d1rVp33DUu2Nw57prP6dAwtmcc38yc4HHbqyqxPSE8vroM0hubvjGBftJ9ENK43dJF2yGBW/jO5AdBO+DXf3h3/Kr3zISM7mMb9DSDurbnq5vGM65XHGE2K4O7Nq4euilGJcVitSiuGJfEe1eP9ri86yGUgh5jzeJoe1d5Pr7SDr88DAn9YNglHg8XbZ/Hrcsao0lbl4nW7acHYPFLMGOhWTGtPktehh/vg14nmDXXG7FzC5ito0JtzTf5xmfPnZ0C702F4hy46BNIGl//sctmmZX3LvoUjjqlcc8nAp6vty4TwnvH3m5mBP70YN33a20GMX+8D/pPhekfNTpoAy0WtJv83LHJcOUPZnG0D86GbT/XfVxpnplpmXQs9JEp4sKQwC18KywGjr8bdv4K23859D6t4cf7TSAadjGc+3Zgrb/S3Np3hSu/N+vQfDQdNnxR+5hFz0FJduN3TBJtkgRu4XujrjFrFv/09/LKploAAAWlSURBVOolPR2V8M3NsPRlOGYGTH3J1IAf6SLi4Yq5Zkbl7KvMRrUueWmw9FUYcoFZTlQIJwncwveCgs1Gxgc3wJqPoKIc5lxjBiOPuxOmPNmw7a3autD2cOnnZuf0r28ywRpg3uPmU8rkB1qydSIASZdH+MfAs8wA5LzHYeNXsO0nOOkxM3lE1BYcAdM/hjlXm82u0zfAmo/N9Yqufz1ycWSSbo/wD9cuPwX7zcDbGc9L0PYkKATOfQeGXgR/vm/GCybc1tKtEgFIetzCf3qMhVP+afLd/U5r6da0DtYgOPNls4BUh/4BsU2WCDwSuIV/jb2hpVvQ+lgs8ulEuCWpEiGEaGUkcAshRCsjgVsIIVoZCdxCCNHKSOAWQohWRgK3EEK0MhK4hRCilZHALYQQrYxfNlJQSmUAuxv58Hgg04fN8SVpW+NI2xpH2tY4rbVtPbTWCd6cxC+BuymUUiu83QWiuUnbGkfa1jjStsY5EtomqRIhhGhlJHALIUQrE4iBe1ZLN8ANaVvjSNsaR9rWOG2+bQGX4xZCCOFeIPa4hRBCuBEwgVspNUUptUUptV0pdU9Lt6cmpdQupdQ6pdRqpdSKAGjPW0qpg0qp9TVui1VK/ayU2ub8GhNAbXtYKbXXef1WK6WafVcFpVQ3pdSvSqlNSqkNSqlbnLe3+HVz07ZAuG6hSqnlSqk1zrY94rw9WSm1zHndPlFKBQdQ295RSqXUuG4tttOyUsqqlPpTKTXX+bNvrpvWusX/AVZgB9ATCAbWAANaul012rcLiG/pdtRoz3HA0cD6Grf9C7jH+f09wP8FUNseBu5o4WvWGTja+X0UsBUYEAjXzU3bAuG6KSDS+b0NWAaMAT4FLnTe/hrw1wBq2zvAuS153Wq08Tbgv8Bc588+uW6B0uMeDWzXWu/UWpcDHwNntnCbApbW+jcg+7CbzwTedX7/LjCtWRvlVE/bWpzWer/WepXz+wJgE9CVALhubtrW4rRR6PzR5vyngcnAbOftLXXd6mtbQFBKJQKnA284f1b46LoFSuDuCqTW+DmNAPnFddLAT0qplUqp61q6MfXoqLXeDyYQAB1auD2Hu0kptdaZSmmRNI6LUioJGI7poQXUdTusbRAA1835cX81cBD4GfPpOFdrXeE8pMX+Xg9vm9badd3+4bxuzymlQlqibcDzwF2Aw/lzHD66boESuFUdtwXMOycwXmt9NHAqcKNS6riWblAr8yrQCxgG7AeeaamGKKUigTnArVrr/JZqR13qaFtAXDetdaXWehiQiPl03L+uw5q3Vc4nPaxtSqlBwL1AP2AUEAvc3dztUkqdARzUWq+seXMdhzbqugVK4E4DutX4ORHY10JtqUVrvc/59SDwBeaXN9CkK6U6Azi/Hmzh9lTRWqc7/8AcwH9ooeunlLJhAuOHWuvPnTcHxHWrq22Bct1ctNa5wHxMHjlaKeXabLzF/15rtG2KM/WktdZlwNu0zHUbD0xVSu3CpH4nY3rgPrlugRK4/wD6OEdcg4ELga9buE0AKKUilFJRru+Bk4H17h/VIr4GLnd+fznwVQu25RCuwOh0Fi1w/Zz5xTeBTVrrZ2vc1eLXrb62Bch1S1BKRTu/DwNOxOTgfwXOdR7WUtetrrZtrvFGrDA55Ga/blrre7XWiVrrJEw8m6e1vhhfXbeWHnWtMfp6GmY0fQdwf0u3p0a7emKqXNYAGwKhbcBHmI/Odsynlasx+bP/AducX2MDqG3vA+uAtZhA2bkF2jUB87F0LbDa+e+0QLhubtoWCNdtCPCnsw3rgb87b+8JLAe2A58BIQHUtnnO67Ye+ABn5UlL/QMmUl1V4pPrJjMnhRCilQmUVIkQQggvSeAWQohWRgK3EEK0MhK4hRCilZHALYQQrYwEbiGEaGUkcAshRCsjgVsIIVqZ/wfWOTVryRfdTwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\zhaon\\Anaconda3\\lib\\site-packages\\numpy\\core\\numeric.py:538: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return array(a, dtype, copy=False, order=order)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXd8W9X5/99H25K87Th2EmJn74QkZBM2JaEQCqGslFFG6ZemFAotbeFHKbSlBcoqZZbVQlhlQxiBDAJJIJPs5SzHdrxtWbIlS7q/P64kL3nEQ9KVz/v1ykuxdK/ukXTv537Oc57zHKEoChKJRCLRPrpoN0AikUgkPYMUdIlEIokTpKBLJBJJnCAFXSKRSOIEKegSiUQSJ0hBl0gkkjhBCrpEIpHECVLQJRKJJE6Qgi6RSCRxgiGSB8vIyFByc3MjeUiJRCLRPBs2bChTFCWzo+0iKui5ubmsX78+koeUSCQSzSOEONSZ7WTIRSKRSOIEKegSiUQSJ0hBl0gkkjghojF0iUTS92hoaKCgoID6+vpoNyXmsVgsDBw4EKPR2KX9paBLJJJepaCggMTERHJzcxFCRLs5MYuiKJSXl1NQUEBeXl6X3kOGXCQSSa9SX19Penq6FPMOEEKQnp7erZ6MFHSJRNLrSDHvHN39njQv6G6vj1fXHWbvMUfEj702v5z3txQil/HrATxO2Pgf2PY/8DVEuzUhFEXhgy2FrM0v7/k3ry2B7/4Nez6DMOfQnmMOlnx7GLfX1/PHlsQlmo+hP7JsL0+u2E+G3cTK20/DZo7MRzpY5uTSZ9aG/j5/Yk5EjhuX1FXB8+dA6U7178FzYNFbYEyIbruA97cUcvNrmwFYftup5GXYeuaNj22Hl84HV5n699Sfwrn/gIBDc7q9XPbMWsqdHo5UuPjNOaN65rh9lMcee4wnn3ySSZMmUVpaSllZGb/73e8oKirihhtuwGq1RruJPYKmHbqiKLy5/ghpNhNltR4+33EsYsd+ec0h9DqB2aDj7Y0FETtuXPLpH6B8L1y6BM7/Jxz6Gj79fbRbBcC7m45iMerQCXh1Xacm63VMQx28dgXojXD9lzDzF7D+edj1YWiTT7cXU+70kGo18sb6I7IX2E3+9a9/8fHHH3PzzTfT0NDA5s2bueSSS3jkkUdwuVzRbl6PoWlB31/qpKzWw+0/GElOsoUPvy+M2LGX7TzGaSMzuXjqQL47UEGDzx+xY8cVZftg8ysw4+cwaj5M/on6//UvQOGmqDatvsHH2vwKLj3pBE4f1Y+Ptxb3zBt//RhUHoAfPQ0DpsCZ90DmKPji3lDo5cPvixiYmsBvzxlFWa2HA2XOnjl2H+TGG28kPz+f+fPnM3v2bDZv3sykSZN49NFHKSws5LTTTuO0006LdjN7BE2HXL49UAHA9Lw0vh/Zjw+/L8TnV9DrencApqi6jsMVLq6cOZiclAT+u/Yw2wtrmDQopVePG5dsehmETnWpQU69A75/HVbcD5e/HrWm7Siqoa7Bx4wh6RRXW1m2s4QjFS4GpXWje15fDd88DqPPhyGnqM/pDTD7V/DujXBgFb7cuXx7oILzJ+UwZXAqABsOVTIk094Dnyq63PPBdnYU1vToe47JSeLu88a2+fpTTz3FJ598wldffcW2bdt48MEH+fBDtTf08MMPs3z5cjIyMnq0TdFC0w59Z1ENdrOBvAwbJ+Wm4qj3sru49wdHgzeSGUPSGZ2dBKgDWJLjxOeFzUtgxDmQ2L/xeUsynHQ97PkESvdErXn7jtUCMKp/IjOHqhd8twdHN7wEHgec/Ovmz4/9EZiTYOub7CyqodbtZXpeGkMz7SSaDWw9Wt2940r6BJp26AfKnAzJtCGEYOrgNAA2Hq5kTE5Srx732wMVJJoNjM5OQlEUTHod+0tqe/WYcUnBt+AsgQk/bv3aSdfB6odh7RNw3qORbxvqTdps0DEozYoAbCY92wtruLirb+jzwrqnIPdkyJnU/DWjBYafBbuXsjVLFfvJJ6Si0wnyMm1xE3Jpz0lLuo+mHXp+aS1DAlkHg9ISSLQY2Fl0fN25fSUOCiqPb1Bke2EN4wYko9cJDHodeRk29klBP372LQOdAYaGiV/aM2H8Qtj6lprS2JtUH4Uj34G7eS9rb0ktQzPt6HUCnU4wJieJbcfplI9UuBp7jQdWQs1RmHZ9+I1HnQuuMpz7v8FuNjAwVc3yGZweP4IeayQmJuJwxE/vWrOCXufxUVhdH4orCiEYk53EjuMQ9Lc2FHDWw6s4/cGVrD9Y0al9fH6F3cWOUKgFYEgcOaiIsvdzGDRdDbGEY9Ll4KmFnR+Gf7271NfAOzfCw2Pg32fCQ6Ng3dOhgcmCSheD0xvj5WNzktlRVIPP37mMkw2HKjjjoZX84JFVvLWhALa+CeZkGP6D8DsMOxOEjuTCrxnVPzE0ySQv3UphVZ3MR+8FbrjhBubNmxc3g6KaFfSjVXWA6syDjMlJYleRo1MXXE19A/d+uIOJA1NIs5n44wfbO5UadrjCRV2Dj1HZiaHnclISKKyuk6llx0NdFRR/D0PauZBOmAUpJ8CWV3v++B4nvHIxfP8GzL5ZTZk8YSYs/Q0s/wuKolBUXU//ZEtol7E5Sbg8Pg6Wd3zzVhSFO9/dTobdxMRBKdz/wSaUne/DmPPU8Eo4LMko/cdzgnNzs/MrN8OGX4EjFXXd/th9lYMHD5KRkcGpp54aGhAFWLx4Mbt27WL58uVRbF3PoVlBL6lR6x1kJTVeHGOyk6hr6NwF997mQqrrGrjn/LHcfOZwth2tYdORqg73C4Z0RvdvdOjZyRbqG/xUuWJnhmPMc3SD+jhoWtvb6HQw8TLIXwnVPZzr/+nv4cg6WPg8nPUnNWXy8jdg0iJY9Xfqtr6Py+Mju5mgqz2J7Z3I0th4uJKdRTXcfOZw7j5vDNM83yI8ThgfZrygCc6saUxU9jCmX6NRyU5W/19cLasVStpHs4J+zBFG0AODoZ1Ji/pgcyEjsuxMGJjMDydkYzHqeG/T0Q7321VUg07A8KzGFLKcFPWCK6yWDqrTFKwHBOSc2P52Ey8FFDVc0VPs+Qw2vAizFsPYCxqf1+ngvEeg/wRMn9xGErX0T24U1mH97Bh0gt3FnTi/thRhMeo4d0IOJw5K4XLrt1Tq0iB3Trv75dsmYBENTDYcCD0XvKkU10hBl7SPdgW9xg00F/Rh/dQBrI5SF51uLxsOV3LWmCyEECRajMwamsGqvWUdHndnsYO8DBsWoz70XPCCK6qSF1ynObpenUxj6SAjKW0IDJiqDo72BF4PfPJbyBgJp9/Z+nW9ERY8gb6unNsMbzZz6CaDjqGZdnYVdTyItu5ABVMHp2E3GxD1VczwbeCdhhk4PO1PQNvgGw5Arnt36Llg2OeYFHRJB2hY0Ouxmw3Ym9RuMRv05GXY2NWBoG88XInPrzA9Lz303NzhGRwoc3Kkov2Ml51FNYzq31yEgg69SDr0zqEoqkMfOLVz20/4MRzbBsd2dP/Y65+Hinw4+z4wmMNvkz2B/YMu4jL9lwxQmpeTGJWd2OH5Ve1qYFdxDdPy1FRadryPQWngHe8svt7Xfh77hgoz5aRiKdsWes5i1JOcYJQhF0mHaFbQS2rc9EtqfUGO7J/I7mPtd4m/PVCBXieYHJiFB3DyiEwAVu0tbXO/WreXgso6RvVPbPZ8us2EEFBW6zmej9B3qciHuorOC/rYH4HQdz/s0lAPXz2k5oEPP6vdTZdnXY0PHVkbH272/Mj+iRytqqOmvu3xkm8PVqAo6gxmALa+iZI2jAOm4XzVzvkFsKvYwdGEEVC0pdnz/ZMsMuQi6RDNCvqxmnqyEltnC4zun8iRijpq3d429113oIJxOUnN3P2QDBvZyRbW7G/bQQVng45sIegGvY6UBCPlTvfxfoy+ybHt6mP2xM5tb+8HQ05Vwy7+btTM2bJEncg09/ZQVcO22F+fyBu6+ei3vgHl+0PPBwfD2wvrfXugHJNBx8RBKVBTCAdXIyZczIknpLLhUGWb+9U3+MgvrcWVPhZKd4GnsbeYlWyhxCHPL0n7aFbQi2uap5QFGdnBBef2+th8pIqTctOaPS+EYFpeGusOVLSZfhh8z5YhF4B0u5kKp3TonaJkJyDUOHZnmfBjqD6szi7tCn4ffPMY5EyGvLkdN9Hh5rPkhWpMfc0/Q88H0wl3tTPfYd2BCiYNSlHHWba9DSgwbiFTB6ex+5iD6rrw7n5fSS1+BQwDTgTFH/ieVNJtJiqkYYg6L774Ir/4xS863KawsLFQ4HXXXceOHT0QLuwEHQq6EOJ5IUSJEGJbk+ceEELsEkJ8L4R4RwgR0apUiqK0GXIJhkPaEvTvC6rxeP2N8c0mTM9Lp9Th5mB5+Dj67mIHVpM+NIOvKcESvpJOULIDUnPBdBxFrkadC4YENW+8K+x8Xw31zPlVh+4coMLpQSRmqWmTm16BWjVU0j/JQpLF0GYc3VHfwLaj1cxoEm4h50TIGMbU3FQURR3DCdvEwE2i37BAWYDSXaHXUq0mKp0yLVYLtBT05557jjFjxkTk2J1x6C8C57R47nNgnKIoE4A9wO96uF3t4nB78fj8ZNhaC/qAlARsJn2bqWXBwlotHTrA9CHqc+vaKMC0q7iGEVmJ6MJUc8ywm6RD7ywlO6HfcZ7g5kQ1V3z7O8e/opGiwOpHIH0YjPphp3apdHlIs5nU1EafB759BlB7cqOyk9oU9A2HKvErMC0vHUp3Q9FmGK9Wf5k0KAW9TrDhYFuC7sBi1DEgbzToTc0EPc1mpNbtlbNFu8EFF1zAlClTGDt2LM88o/6edrudP/zhD0ycOJEZM2Zw7Jg6CP7BBx8wffp0TjzxRM4888zQ80EcDgd5eXk0NKjnYk1NDbm5ubz55pusX7+eK664gkmTJlFXV8epp57K+vXrAfjkk0+YPHkyEydO5Iwzzujxz9hhcS5FUVYJIXJbPPdZkz/XAgt7tlntUxVwKilWY6vXdDrByP5tZyKsO1DByKxEUm2mVq8NybCRYTez7kAFl047odlriqJO+f/B2P6t9gPVoZfXyi5xh3jdUL4PRp93/PuOv1hdom7/lzCijenz4chfoQrreY+BTt/h5qA69FSrCTKGq72D755V3b3Jxuj+ifxv41EURWm1BuS6AxUYdILJg1Ng1RPqYG5A0G1mA2Oyk1h/KHyZiV3FNYzMSkRvMEL6cChrrDSZFjAvlc4G+id37jPEJEvvgOKtPfue/cfDvPs73Oz5558nLS2Nuro6TjrpJC666CKcTiczZszgz3/+M7/5zW949tlnufPOO5kzZw5r165FCMFzzz3H3//+dx566KHQeyUmJnLqqafy0UcfccEFF/Daa69x0UUXcfHFF/PEE0/w4IMPMnVq80H/0tJSrr/+elatWkVeXh4VFZ0rN3I89EQM/afA0rZeFELcIIRYL4RYX1ra/gh/Z6l0qU441dpalEGNo+8qdrSKhXt9fjYcrAgbbgm0lel5aSEX35RSh5tKV0OrAdEg6TYzVXUNeOVCF+1TthcUH/Qbffz7Dj0DElKPP9tl9cNg7x+YpNQxDT4/jnqv6tABZv0S6iph038B9fwKZjy1ZF1+ORMGJmM16NSa7sPOUAd1A0wZnMrmI1V4vM3PE0VR2FlU01gjKHNEK4cOyF5gN3jsscdCTvzIkSPs3bsXk8nED3+o9tqmTJnCwYMHASgoKOAHP/gB48eP54EHHmD79u2t3u+6667jhRdeAOCFF17gmmuuaff4a9euZe7cueTl5QGQlhZeh7pDt8rnCiH+AHiBV9raRlGUZ4BnAKZOndojxU5Cgm5r7dBBjaMv+fYwx2rczQZOdxTV4PT42hR0gGl5aXy0tajVQgZBx9+moNtNKApUuhrITGwjv1nSKFJdEXSDSU1h3PIauGvB3IkFHwo2qFUOz7q37bzzFjSeXwFBP2E6DJqhLkwx9aeNA6PFjmbniMvj5fuCaq47eQgc/EqtrHj2vc3ee1peGi9+c5BthdVMPqExbTZoGEIpsZmjYPu76nJ1xoSQedG8oHfCSfcGK1asYNmyZaxZswar1cqpp55KfX09RqMx1MvS6/V4vWp23OLFi7n11ls5//zzWbFiBX/84x9bvefs2bM5ePAgK1euxOfzMW7cuHbbEK5H19N02aELIa4CfghcoUS4KlWwZkpKGw49eFHsahFHDzrv9gQ9FEdv4dK3FaplU0eHyXAB1aEDMnWxI8r2AALShnZt//EXQ4MLdn/cue1X/wMsKTC1fffUlODgY1rT8+vkX0P1Efj+dUZmhc902XS4Cq9fUc+hjS+plRVHzm+2TfDcW5ff/PwKVgkdFXToGSMARe3RoBoGgAqXxgU9SlRXV5OamorVamXXrl2sXbu2w+0HDBgAwEsvvdTmdldeeSWXXXZZM3feVknemTNnsnLlSg4cUMs6xEzIRQhxDvBb4HxFUSK+wmpHIZfgRdGydvXa/HIGp1ublQtoyYh+iaRYjXx7oPnA6KbDVeRl2MLG3oFQ97xCZrq0T0U+JA9qu+JgRwyaAal58N1zHW9bsktdeHn6z9RB1c420RmmBzj8LDVv/quHsBnghDRrq3Gadfnl6ASclFYPO96DExeBsXlGVIbdzLB+dta1OL82H6lCJxrrEZE5Sn0sVUsAhBy6HKfpEueccw5er5cJEyZw1113MWPGjHa3/+Mf/8jFF1/MySef3O7ydFdccQWVlZVcdtlloeeuvvpqbrzxxtCgaJDMzEyeeeYZLrzwQiZOnMgll1zS/Q/Wgg5DLkKIJcCpQIYQogC4GzWrxQx8HuhCrFUU5cYeb10bVDo9CAHJCeFDLskJRkZmJfJtk2wCr8/P2nx1ncb20OkEJ+WmNXPoiqKw6XAVc0e0/cMGHVSZ1rvEvU1FPqTldX1/nU5dRHrpb9RFKQad1Pa2K/8GRitM+9lxHSJoGNKa3ryFUCckvb4Itr/NhIHDWH+wslk3evW+MsYPTMG+9T9q3vu068K+/7S8NN7f3Hz92/UHKxnZP4kkS+CcTh8KCKhQJzWlWNXZyBWyomeXMJvNLF3aeqivtrZxYZqFCxeycKGa37FgwQIWLFjQavurr76aq6++OvT36tWrWbhwISkpjZnbF110ERdddFHo7xUrVoT+P2/ePObNm9edj9IuHTp0RVEuUxQlW1EUo6IoAxVF+beiKMMURRmkKMqkwL+IiTmoceoki7HdxaBPyktl46HKUG30LQXV1Lq9zB7a8WKw0/PSOFTuCtXOKKiso6zW3Szm2ZKUwM2lrUkjkgAV+QGx6gaTrlAXxfimnaXpCtbD9rdh5k1gS297u3BNDNyU01r2AEeeC1nj4Mt7mT3YRnFNfWjOQnVdA5uPVHFWnkVNcRxxjlpYLAzT89KodXvZHgjjeX1+Nh2uZGqTUhQYzJA0ACrU7rleJ0hJMMrJRTHE4sWLueOOO7jrrrui3ZQQmpwpWunykBomZbEpJ+WqF01wssaqPaUIATOHdnxxB4t2BbvFawJ56VMGty3oSQFBr5GC3jauCjVbpA2h6zRmO0z/Oez8AI6EmTnq98End4AtU1284jipDAh6qzEanQ7O+StUHebsanWCU7BUxJr9ZfgVuLD+baivglPvaPP95wzLQCfgs+1qbvPOIgdOj4+puS3Or7Q8qGwso5tqk3MdYonHH3+cffv2MWLEiGg3JYQmBb3K1dDmgGiQmUPTEQKW7TyGoih8tLWIkwanNe9Gt8GYnCSSLAaW7yoBYNmOY+QkW1oV5WqKxajHbNBJQW+PgNvstqCDOuHHngUf3qrmtjflm8eh4Du1ouJxxM5DzXR5sJsNmAxhLo+8uTD2R6RteIw59qOs2qOm4n68tZjJlmKytz0N4xa2XgS6Cel2MzOHpvPx1iIUReGzHcXoBMwe1qL3mJrb+J2h9gJr6tquURTLyNW8Okd3vydNCnpnHHq/RAsz8tJ5Z9NR1uZXsK+ktsP4eRC9TnD+pByWbitmZ1ENK3aXcvbY/h2mHCUlGGXIpT0q8tXHrma4NMVsVycKHdsK7y8GX0Dotr4Fy/4IYxbAhK4NOlXXNbQ5PgPA/IcQ1gyeEA+wf9dm9hxzsGH7Lp41/wNhSYJ5f+vwGPPHZ5Nf5mRtfgXvbS5kel46GfYWaZVpeWoxMbca501KMLZb5TFWsVgslJeXS1HvAEVRKC8vx2LpYsIA3cxDjxZVroZQ6lh7XDUrlxv/u4Ern19HZqKZiyYP7PQxrp6Vy+vfHeH8f67G51e4ZnZuh/ska/SCixgV+wGhOs+eYOQ56iIVX96nzj60ZcCBVWomzAVPdapmSzgc9V4SLe1cGrZ0uOJN7C/+kPcNd7D6n6/wgX43KT4fLHpfbUcHXDBpAA9/vpcrn19Hg0/h//0wTCmE1MDgceVB6D+OJIuRgxpcjHzgwIEUFBTQUxML4xmLxcLAgZ3XqZZoVNA9JHfg0AF+MDaLn87OY9XeUu4+bwwJps5PmR7WL5E/LRjHM6vy+ensXAan2zrcJ1k69PapOKAO9HU1ZTEcc29XHf+6p8BZDqf+Hubcok5C6iI1dQ2hMZE26T8O/Y1fUfTm7xlRuBFnv2mkXXgfZI3t1DFsZgP/+PFE/vThDuYOz+SM0f1abxTMBqo8oAp6goGaeu2FXIxGY2h2pKR30Zyg+/wKTo+PREvHgi6E4P+d1/UqZ5dNO4HLWtR0aY8ki4FSmSfcNlWHIaXz32enGXeh+q+HcNR7yUnpxE0nZRBDrv9Pl48zd0Qmy249pe0Ngg49EEdPsqiGIRIzDiXaRHMx9NqAQ0lqr0scJZI1PGgVEaqPQMqgaLeiQxzuhk4Zhl4nIUWtXRPIdElOMOLzK7g8suKiJDyaE/RgjLrdGGeUkCGXdvB51dV7kmNf0GvqvLFjGFIGQ9URoElqrBynkbSB5gTdEXDoMeGgWhDMQvD75Wh+KxyFapXFGHfoiqLgqI8Rhw6QPBCqCwBCs0hlL1DSFhoU9Nh26IoCtR55wbUi4DJj3aE7PT78SgydX8mD1FCVopCUoLZJOnRJW2hO0IOLP8eMg2pCsEtcLetttKbqsPrYG4OiPUjQMHSY5RIpkgeCpxbqq5o4dHl+ScKjOUFvDLnEiINqQuiCkw6qNdVBh971HNtIEHPnV/D7qi6QMXRJh2hQ0GM75AKyQFdYqg6DrV+rcrKxRtD9JsVKDzAYoqouCA3Uyh6gpC00J+g1obTFGLngmhC8ydRqcPJHr6OVlMVYc+gpjYKeGOoByvNLEh7NCbqj3otRLzCHK5wUZezmgKC75QXXiqojMT8gCk3TYmPEMFgzQG+G6iOYDDoSjHoZQ5e0SeypYgcEU8picaacPeDqnFLQm+P3q6l3GnDooR5gQow4dJ0Okgc0pi4mGGQMXdImGhT0DgonRZGgQ3dIQW+OsxR8bkiO7QwXaJLlEisOHdSB0UDap81swOmWM0Ul4dGgoDfErKCbDToMOiFj6C2pUd0lyQOi245OEJMhveRBIYeeaDZIwyBpkxg6azuHo95LojmG3FMThBDYLQYZQ29JTZH6mNS5evTRpKZOXd4wpkJ6yQPBUQS+hoBDl+eXJDzaFPQYdeighl2koLfAERD0xNgX9Jg8v5IHAgrUHMUuBV3SDhoU9IbQ4GMsYjcbZMilJTWFoDOoa3zGODFVxyVIsGdTU4TdbAilVkokLelQ0IUQzwshSoQQ25o8lyaE+FwIsTfw2PbqyT2Mo94bWwNWLZAOPQyOIrD3VzM2Ypyaem/sZLgEScxWHx1FashF1gqStEFnrrAXgXNaPHcH8IWiKMOBLwJ/9zp+v0KtJwa7xE2wW2SXuBU1hZCUHe1WdIraem8oWylmaCLodovaA5Trc0rC0aGgK4qyCqho8fQC4KXA/18CLujhdoXF6fGixFIlvDDYZRZCaxxFjaIU4zg9XmyxJugJqerkIocacvH6Fdxef7RbJYlButoHzlIUpQgg8BhmQUQVIcQNQoj1Qoj13V0kNpZroQeRMfQw1BRpIsMFwOXxYTPFmKALofZwAjF0kJPXJOHp9aCmoijPKIoyVVGUqZmZ3RsUC8amY65L3ASZhdACtwM8Du04dLcXq7nzi4lHjMRscBSHeg9ynEYSjq4K+jEhRDZA4LGk55rUNk4tCLrFgNPjwydXLVLRUA661+fH7fXHnkOHgKAXynpBknbpqqC/D1wV+P9VwHs905z2CS6OazXFoIMKEOoSy0wEFUeh+qgBh+4MnF8xF0OHkEO3B859GdaThKMzaYtLgDXASCFEgRDiWuB+4CwhxF7grMDfvU7QocfkBRcg5KDkBaeiIYfuCtyEbbFoGJKyocFFos4FSMMgCU+HyqgoymVtvHRGD7elQzTh0GXFxeZoyaEHil5ZY9EwBL6/FG8ZALWyQJckDLE/06MJQVeiBYcuUxcD1BSBJRlM1mi3pENi2qEHBN3uCQi67AFKwqAtQZchF+3hKNJEDRdo4tBjclC0PwBWt5p/IHuAknBoTNDVCy7BGIMOKoAMubRAQ7NEGw1DDJ5fAYdurlMFXfYAJeHQlKC7PF4SjHr0uhgqbdoCGXJpgZYceiyH9ExWsCSjcxRhNemlYZCERVOC7vT4YtM9NUGGXJrg90HtsVC4INYJDrrHZB46qDfGwPR/eX5JwqEpQXe5vbEZ32yCnMnXBFc5KH7NCHrQ9cbkTFFQv8egoMu0RUkYNCXoTo8vplMWAYx6HSaDLuT2+jS1x9RHe1Z029FJQmmxsTpGk5SjTi6SFT0lbaAtQXfHYGnTMNhM+lAKXJ9GY4Lu9HjVdWH1MXpZJPZXBd2okyEXSVhi9MwNj9Pji81JHy2wmuTK7ADUBkr82NssxhlTON0xWDq3Kfb+oPjoZ3TJkJ4kLJoSdJfbG5uTPlpgM0uHDjRx6NoQdJc7xgfdA99jtq5aTv2XhEVbgu7xxfygKAQcuoyhg+MYmBLBZIt2SzqF0+ON3QwXCIWuMnXVuGQPUBIGTQm6uppMDDuoADazHpfsEqsOXSPuHIKGIYbPr8B3mUGVDLlIwqIpQXe5pUPXFLUlmhkQBS3E0NXvMtVfidvrx+uf4AuGAAAgAElEQVSTy9BJmqMZQfd4/Xh8fuwacOhWk546GeOUDr2nMdvBaCPFry7xK02DpCWaEfTgIKN06BpCYw691h3jMXQAez8SvQFBl2EXSQs0I+iNq8nEsIMKYDPJGDoNdeCu1pxDj+mQC4A9C1uDKugyk0rSEs0IelAgNeHQzQZcDT78fXld0VAOunYceswuEN0Uez8S3HKRC0l4NCPoWnPoigL13j58wWlM0GN6geim2LMw16uC3ud7gZJWaEbQtebQgb49W1Rrk4oaYn95QwDsWRg9VZhokKmLklZoRtCDJ69WarlAH49xaqyOS3CiTuzH0NUbZDo1sgCcpBWaEXQtLBAdJNiL6NsOvQQQYMuIdks6RW2oBxjj51dwtqiQk4skremWoAshbhFCbBdCbBNCLBFCWHqqYS2J6dVkWhCM8/d5h25NB70x2i3pFMHfKuZ7gAGHnimq+vb5JQlLlwVdCDEA+CUwVVGUcYAeuLSnGtYSl1uDDr0vd4k1loMe0wtENyXk0KtlloukFd0NuRiABCGEAbAChd1vUnicGppYFHLofblLXHsMErUj6K5QDzDGDYMtE4AcfU3fPr8kYemyoCuKchR4EDgMFAHViqJ81nI7IcQNQoj1Qoj1paWlXW6oy+PDYtTF9ALRQWzSoWvPoXs04tANJkhII9tQI0voSlrRnZBLKrAAyANyAJsQYlHL7RRFeUZRlKmKokzNzMzsckNrNbJaEUBCICzUZ+u5KIr26ri4NeLQAexZ9BNVMuQiaUV3Qi5nAgcURSlVFKUBeBuY1TPNao0WFogO0ucden01+NzSofcW9n5kimoZcpG0ojuCfhiYIYSwCiEEcAaws2ea1RotLBAdxGLUIUQfjqFrbJYoNBa60sKKWNizSFNk2qKkNd2Joa8D3gI2AlsD7/VMD7WrFS5PjNeqboIQAltfrriosVmioIEFopti70eqv7LvGgZJm3RLIRVFuRu4u4fa0i5Ot49EizYEHdT0yj6bJ6yxWaIQXE9UI+eXPQuzUo/f7Yh2SyQxhgbsiIor1td7bIHNbOi7M0VDIRdtOXSthPSCN8oET1mUGyKJNTQj6E4tOSikQ0dvAktKtFvSaVxun3YMQ+BGafVURLkhklhDO4KukQWig9hMfdyh27NAxP6cgSBOjwZqoQcJOPREbzmK0odr7ktaoRlB18oC0UGs5r7s0Is1FW4BNctFK/McgoKeQTV1DX3UNEjCoglBDy4QrYmUsgB9O8tFW7NEQQMLRDclIRW/MMiKi5JWaELQ64KTPrTioAjE0PvqxaaxWaIQCOlppQeo0+E2p5NJdahonUQCGhH0UOlcrTgoAlkufdGh+7zgLNOeQ3f7tBNDBzwJGfQTldKhS5qhDUF3a6cWepAEkz7Us+hTuMoARTr0XsaX0I8MUS1XLZI0QxuCrqEFooPYTHo8Pj8erz/aTYksGpxU5PMr1Df4NTXortgyyRTVIbMjkYBGBF1LC0QHCba1z7l0LdZx0Uot9KYkZpFBNU63J9otkcQQmhD0kEPXkKAHxaHP1awOOnRb10slRxrNLBDdBH1SFgbhx+uQs0UljWhC0IP53FoatAo69D6Xi67BkEvjaljaOb+MSf0BUII9IokEjQh6cMalJh16X0srqy0FcxKYrNFuSadxafD8MqVkAyCkoEuaoBFB116Ms3Gh6D7o0DUUboEmDl1D55cxWRV0XZ0UdEkj2hB0DS0QHSTo9vrcxA9NzhINznPQzvkVTAs11ckYuqQRTQi6lhaIDmLty4OiGstBr9XgoCgmO3WYMddLQY91Kp0e3tt8lGM19b1+LE0IutOtrUkf0MSh97W0RacGHboGQ3oIQaUuFausiR7z5JfVcvNrm9lV3PsLkmhCJW89awTXzsmLdjOOi5BD70sTPxrq1QWi7VqLoWtogegm1OhSsTXImuixTmNSR+8bBk2cwel2M+l2c7SbcVxYjeqP16cmFjm1N6kImk5c05BDBxzGNLI8BdFuhqQDIlm6RBMhFy1i0OswGXR9q0BXban6qDFBd3p8mAw6jFpYILoJLmM6yX7p0GOdSE6M1NYZrDFsfW0ZutCkIm0Niqrr1WrLnQPUmdNJVhzgldP/Y5lITozslqALIVKEEG8JIXYJIXYKIWb2VMPiAavJ0LfKm4am/WtL0GvdXs3FzwHclsBYhbM0ug2RtEskJ0Z216E/CnyiKMooYCKws/tNih/sZkPfykMPzlrU2MQil9unneXnmtCQEBR0ObkolnF5vAgBFmPvB0S6fAQhRBIwF/g3gKIoHkVRqnqqYfGA1azvW3nozhJISAODKdotOS40tUB0E/xWVdAbqouj3BJJezjdPmwmAyICi6Z355YxBCgFXhBCbBJCPCeEsLXcSAhxgxBivRBifWlp3+oa2kyGvpWHXntMcwOioM4V0No8ByD0XXuqiqLcEEl7ON3eiGVQdUfQDcBk4ElFUU4EnMAdLTdSFOUZRVGmKooyNTNTW13x7mIz6/tWHnptieZy0CGyF1xPogsIurdGOvRYxunxRiyk1x1BLwAKFEVZF/j7LVSBlwSwmQx9K+SiUYfu9Hi1Ne0/QILVSpViw+84Fu2mSNrB5YncerVdFnRFUYqBI0KIkYGnzgB29Eir4gSrWd/HBkVLNSnoLrdPW9P+A1jNesqU5MbBaElM4oxgFlV3j7IYeEUIYQLygWu636T4wWbuQ2mL7lpocGouwwW0t0B0ELvZQKmSQqbMcolpXB4fGfbIJAp06yxWFGUzMLWH2hJ32EwG3F4/Xp8fg8ZmIR43GlypCMDr82tugeggVpOe/SRjqCuMdlMk7eD0eDnBHJkFX+JcZaJLcKCtT0z/Dy0Ora1JRa6GYOlc7YVcgg7dVNe3sse0hsvti9hMZCnovUhwZLtPTP/XbGEuDdZCD2A1GShVkjH6XOBxRrs5kjZwuiM36C4FvRexBn7EPrGuqEYdem0EK+H1NEGHDsiB0RhFUZSIjtFIQe9Fgt2sPpGLXnsMhA6s6dFuyXHRuPyc9kIuFqOOMpLVP6SgxyRurx+/Ern1aqWg9yJB19cnctGDi0PrtCWMtW7trVcbRAhBrTFwA62VueixSKgWunTo2qdPLRRdW6q5cAs0/jZaLM4FUGfKUP8jBT0mcYVWw5IOXfP0qYWia49prmwuNP42WizOBeAxp+BHJ0MuMUrw/JKDonGAva8NimoswwUafxutOnSrxYxDlywdeowSPL+kQ48Dgj9i3KctKoqatqjFkItHm+uJBrGa9FTqUqVDj1GCMXQtFOeSdEBwoC3up//XV4HPo0mHruVBUVCFoowU6dBjlEbDIAVd8+h1ggSjPv5roms0Bx3UQasEox69rvcXH+gNbGYDpcgCXbFKaPk5mbYYH/SJmugaXRwagrP4tBluAdX5HfMlqyEvRYl2cyQtkA49zrCaDH1A0LU57R8iW9q0N7Cb9RT5ktSQV71cATLWCNZxkg49TrCZDfFfnEvLDt3j0+S0/yBWk4FCb5L6hwy7xBwud2CBaIMU9LjAZtLHf5aLowgMFrCkRLslx43T7dXktP8gdrOBUoL1XOTAaKxR6/ZhNerRRWiMRgp6L2M1G6iN9zx0RzEk9ocIrGre02jeoZv1lCqynkus4orw8oZS0HsZu1mPK95j6I5iSMyOdiu6hEvjg6JqxcWgoEuHHmtE2jBIQe9lrCZD/KctOopUh65BtD4oajUZqMGGX2eSgh6DuNzeiE5ak4Ley9hM+vifWKRhh+70+DQ77R+C2ROChoQMGXKJQSK9Xq0U9F7GZjbE96Co2wGeWk2mLCqKEnDo2g65ANSbM9QbqySmcHl8ES38JgW9l7GZDTT4FDxef7Sb0js4At18DTp0j8+P169oe1A04P6c5n5q6EsSU6hZVBpy6EIIvRBikxDiw55oULxhjfdVi4IiosEYemg90Thw6A5TJtRIQY81nG6f5mLoNwM7e+B94pK4X7Uo2M3XoEMPFebSskMPdOerDRngrpaLRccYTi2lLQohBgLnAs/1THPij9CqRfGa6aJlh+7Rdi10aDy/KnWBpehkHD1mUBQFl8cX0bTY7jr0R4DfAG0GiIUQNwgh1gsh1peWlnbzcNoj6KDiNtPFUQxGK5gTo92S46axdK52Qy56ncBi1FEeFPSawug2SBLC7fXj8ysRTYvtsqALIX4IlCiKsqG97RRFeUZRlKmKokzNzMzs6uE0S9D9xe26orXanSUazD7SskMHtf3HSFP/kAOjMUNthBe3AOjOkWYD5wsh5gMWIEkI8V9FURb1TNPig9CgaDzH0DUYP4emy4NpW9BtZgPFfqv6h3ToMUNtfeQFvcsOXVGU3ymKMlBRlFzgUuBLKeatCcY44zrLRYPxc2j8TbQ89R/UG1K51wwmu4yhxxBBh55o0YCgSzpHY5ZLHIZcFEXTDj0aXeLewG4OVPRMzAaHdOixgiPo0CMo6D1yJEVRVgAreuK94o2gWAS7X3GFuwYaXJp16CFBj+AF1xtYTQaqXB5Iypa56DFEyKGbjRE7pnTovYzFqEOvE9S6G6LdlJ4n2L23a1PQa+obMBl0mCO0+EBvYTcbVPFIzJGDojGEo1695mXIJY4QQpBoMYS6X3FFcAAuSZshF0e9lySNu3NQB95dHp/6OziKwB+nZSY0RjR6gFLQI0D8CvpR9TFpQHTb0UVq672aj5+DOk6jOvRs8HvBVR7tJkloEkPXQpaLpPMkmo2h7ldcUR0U9JzotqOLOOobSLRELr7ZW9jMqkNXgmMZcmA0Jqh1ezHqBWZD5GRWCnoEsMetQy8AWz8wmKPdki7hqPdGNL7ZW9jMBnx+BY81IOhyYDQmCBoGEcFJd1LQI0BSvAp69VFI1ma4BVQHFQ8hl9BsZHM/9Qnp0GOCaIT0pKBHgESLEUc8ZrnUHNVs/ByCDl37IZfgTFeHPg0Q0qHHCNEwDFLQI0DcDopWa13QG+Ii5GIPzHR1+oQ6MBocrJZEFUe9N+JzHKSgR4CgoCuKEu2m9Bz11eBxaDbkoigKte74iKFbm5aXSB4IVYej3CIJqA490mmxUtAjgN1sxOdXqGuIo+n/oRx0bQq60+PDr0R20kdvESwvUev2QsogqC6IcoskEHDoMuQSfwRFI66m/wdTFpMHRrcdXaSxEp72Y+jB4mIuj0/9PWqOyslFMUCtW4Zc4pKgoNfEk6DXBFygRh16NKZl9xbBip61bi8kDwKfB5wlUW6VRM1yiaxhkIIeAZICmRRxNbmo+igInWYrLQZvrvEg6I2LqAQEHaDqSBRbJHF7fXh8/oifX1LQI0DwR42rTJeao2pRLr02BTEatap7i+Ayh06PT42hA1RLQY8mjigZBinoEcAej4JeXaDZDBdoGnLRfgzdbNBj1IvGLBeQgh5lorFaEUhBjwhB0YirErpVhyBlcLRb0WWidcH1FlaTQRV0SzKYk2WmS5Rp7AHKGHrcEXchF59XjdGmalfQo9Ul7i3sZkPjqljJA2UMPcrUBHqAkV7eUAp6BLCbDAgRR1kuNUdB8WnaoTvqGxCiMUNE61hN+sZ1a2UuetSpqVN/i+QE6dDjDp1OYDcZ4ifLpfKg+piaG81WdAuH24vdZECni1wlvN4kKcEYcoUkD5Qx9ChTXecBIMVqiuhxpaBHiLgqoVt1SH3UcMilpi4+pv0HSUkwUuUKCvogqK+C+proNqoPU12n/hbSoccpSRYjNXVx5NCFHpK0OUsUVAcVaffUmyRbjSERIS1Pfaw8EL0G9XGqXA0YdAKbSSMxdCHEICHEciHETiHEdiHEzT3ZsHgj2WqkKm4E/ZDarddoDjpApauBVJv2UxaDpCSYqA469LQh6mNFfvQa1MepqmsgOSGyi1tA9xy6F/i1oiijgRnATUKIMT3TrPgj1WqkyuWJdjN6hsqDmg63AFS5PKQkxJFDTzDicHtp8PkhNeDQpaBHjeq6BpKtkTcMXRZ0RVGKFEXZGPi/A9gJaHemSS+TajU1xji1TtUhTQ+IQvQuuN4iJfBZauoawGwHe5YU9ChS7WogJcLxc+ihGLoQIhc4EVgX5rUbhBDrhRDrS0tLe+JwmiTZqg5aab4mursWnKWaTllUFIWqKF1wvUVQ0ENhvbQhUCFj6NGiOhByiTTdFnQhhB34H/ArRVFaDasrivKMoihTFUWZmpmZ2d3DaZZUqwmPz6/9mujBxRM07NBr3V68fiUkgvFAUDyqmsbRpUOPGlVRGnTvlqALIYyoYv6Koihv90yT4pOgG6zUetglKBLBOK0GCYpePGW5BD9LMP+ZtDxwFIHHGcVW9V2qXBpz6EIdvv03sFNRlH/0XJPik+AFV+nU+MBo2R71MWNYdNvRDYLpfXEVcgnn0KFxEpgkYvj8Co56r7YEHZgN/AQ4XQixOfBvfg+1K+4Idu+rtZ66WL5PLZtrSY52S7pMfDr0NgRdhl0iTk2UJhUBdDmRWFGU1UB8zJuOAKlBh6711MWyPZAxPNqt6BbB3yCeYujBqn6hQVGZuhg1Qj1ALaUtSo6P1JYOSosoSkDQR0S7Jd2iKooXXG+h1wmSLAaqg4YhIQWsGWqPShJRygNh1VStDYpKOk9ySNA759BLHe6uxdsVBcr3w9GN4HYc//7t4SyF+uqYdOhVLg8ljvpObRsUvWh0iXuTNJspJCYA9BsNJbvCblvh9FBc3bnvqxU1RVCwAZzlXdtfgzjqG8gvre3UtuW1bgAy7ObebFJYtDt3W2OYDXqsJn2nslz+s+Ygd7+/HaNexwMXT+T8iTmdO8jeZfDZH6A0cBEbLHDSdXD6XWC0dL3xQUIDorEl6O9tPsrtb36PT1G45/yxLJrRfo58Wa0Hu9mA2RDZOhu9TWaimbKAmACqoG9eot7km0xB/2pvKT/7zwZcHh83nzGcW87qZI+r4gB8fBvsW9b43Ih5cM5fG+vHxCGbDldy7UvrqXB6uHLmYP60YFy725fVqjfVjETp0OOadLspdPdui/zSWu75YAczhqQzJieJ297cwpEKV/tvrCjw+d3wykWg+GH+g3DJf2HcRbDmn/Dy+eqEoO5Stld9jKGQS35pLb9563vGDUjipNxU7vlge4dOqrTWTb/EyLun3ibDbqbU0eT8yhwFHkez2ui1bi+/fmMLA1ISmD++P49+sZflu0o6fvPCzfDsaXDkWzjtTrh0Ccy9HQ59DU+dDLuX9sInij71DT7+75WNWE16Lpo8kJfXHGLZjmPt7hO8qabbIn+OSUGPIJl2MyWO9gX9udUHMOgFj156Ik8tmoJOwN8/3d32Dn4/fPBL+PoRmPpTuHE1TLseRp8HF/wLFr4ABevhzavVbbtDyQ4w2WOqyuKjX+xFJwRPLZrC45dNRgjBs1+1P0OyzOEmIw4FXXXoTUMugdJKpY1hl/9tKKDE4eavF47n4UsmMTTTxp8/3onf384MZkcxvHIxmBLhZyvhlNth1Hw4/U74+ddqCutrV8C2+JuK8s6moxRV13P/hRO4/6Lx5KZb+deK9sclymrdJCcYMRkiL6/xH3JRFPA1gN7YrNvZa1Tkw7b/wYFVUH0UDGbVKY2cR459ELvL2hb0+gYfH2wpZN64bDIDgnPtnDyeWL6fW88aQV6GrfVOy++DjS/Dyb+G0+/iu0OVPPbFZmrqGrjgxAFcNfNH6Fzlald53ZMw86YOP4KiKOGrxBVvhaxxoGv7RG1z33bw+xVeWnOQdzcdJdlq4uYzhjFlcFqH+x0oc/LBlkKuP3kI/ZLUkNKCiTm8v/kod583BosxfEiltNbN6Oyk42qjFsiwm6mua8Dt9anhpH6j1BdKdsDwswB4/bsjTBiYzNRc9fu9+cwR/HLJJj7bUcw547Jbv6nfD29dC55auOp9DvizePDVjRwsczJ3RCaLTx+G9aoP4JUfw9s3gDUdhpzSbjuP+xxxFMOW19RrquoQ6Axq2G/YWTD2gk6l0C7dWsQL3xxEURSumpXLDyd0HMb0+xWe/SqfcQOSmD00DeH3smjGYO77aCd7jjkYkZUYdr+yWjfp9uikxMafQ1cUOLwOPrwVnpgB92bCfZnw10Hw7Omw/K/qoGFP43HC5/8P/nkSfHmfOniYPUGdIn94Dbx9PX85sogTa75U2xiG5btKcNR7uXByY42zq2blYtQLXvrmYOsdNrwEXz0EU66G0+9ixZ5SLn92LftKakEI7vlgB7e+sRn/lGth5Lmw7J5263u8s6mAM/+xkmF/WMrCJ79hw6GKxhf9fijeBv3Hh913w6EKLvzX1wz9/cec/uAK3tt8tBNfmnrR3PLGZu75YAcIwZ5iB5c8vZZVezqu+/PMqv0Y9TquO3lI6LlzJ2Tj9Pj4Zn9Zm/uVOtxkRmHAqrcJDsKVB116QiokZkPJTgCOVLjYUVTDeU3E7Nzx2QxOt/LkyvzwdYY2vwKHVsO8v3FQdwI/+tfXrNpTSpLFyFMr97PouXU4FAtctgTSh8Hri9TzJAzr8stZ8M/VDPvDUs5+eCUffV/U/gfyOOGLP8EjE2DZ3RQczufj0gxWlCVRc2Cj2jN9eDys+Bt42g5LPrF8Hz9/ZSPltW6qXA384tVNPPdVB+mcisL6VR9xXeUjLHEvRtyn6shPV5/C/0x3U/HxfW1O2ipzeKIyIArxJOiKArs+gidnw/Nnq3f05IGqIz3tTjjxCvXOvurv8PgUtYt4bEfPHHv3UnhiOnz9KEy8FG7ZDj9bBRe/qJ7ot+yARf+j3pTO33kE32uLVMFvwaq9ZSSaDcwckh56rl+ihfMm5PDm+iONS4wB7P8SPrwFhp0J8x+iuMbN4lc3MbxfIp/8ai7v/t8sbj1rBO9uLuThL/bCuQ+pn/+zO8N+hL8u3cktr28hwajn2jl5HK2q45Kn1/Lm+sBSZlUH1XhsGEFfurWIHz+9luLqeq4/eQhWs56bX9vMXz7e2WExsieW7+O9zYVqW/9vFp/eMpfhWYnc9MpGiqrr2tyv0unhnU1HuXDygFBvBmDW0AwSzQaWbi0Ou199gw9HvbfZPvFC8DO1Ghg9pgrsZ4HY71ljskIv63WCa+fkseVIFRsPVzZ/w7pKWHY3DJqBf+IifvO/7/H5FT5cPIclN8zgySsm831BNbe8vhm/ORkWvQUmG7z6YzUTpgnvbjrKZc+updLVwHVz8tAJwU2vbuT+pbvCnyPl++HZM+Crh6geci4LxKPM8/yN7076By+e8GcmVP2dP+f8E9/gObDiL/DkLNi/vNXbfLKtiAc+3c0Fk3L45FdzWXrzycwf35+/fLyTNfvDZOmEdGQW01ZcwQLDWmzZI2DWYjjtTnQTLyXRJJh26Gl47ERVR1rcwMpqo2cY4kPQD66Gf58Fr10OPjec/zjctkc9wc66R435zfsbXPsZ3LoTTr4VDnwFT82Gj37d9fSr6gL1B11yqXoiX7MUFjyh3kiaotPBsDNZccrr/KXhMnR7lsLTp8Cx7c02W5tfzrS8NAz65j/L1bNzcXp8vPFdQFyLt8HrV6oX68IXUHR67nx3Kw1+P08umhwqrL/49GFcMnUQj3+5j6WHUL+HXR+q31cTnl2Vz9Mr81k04wTevWk2v58/ms9umcvMoenc/tb3vPbtYTXcAmqvowmf7zjG4iWbOHFQCp/eMpffzR/NezfN4SczBvPMqnyeXNl2b2j57hL+sWwPCyblsPj0YQghSE4w8tSiyXj9Cr9/e2ubN4TXvjtCfYOfq2c1z64wGXScMjKTVXtLw+4bHDSMT4eudvObDYzmnKgal4Y6VuwuYXg/O7ktQncLpwwkOcHIcy3HHlY9qIr6uQ/yyndH+PZABXedO4bB6er+54zL5s5zR7NsZwmPf7lPPe8vf0M1K69eHEqbXbmnlNve3MK0vDSW3nwyv5s/mg8Xz+Hy6Sfw1Mr9PPrF3ubH3bdMHYCtLabgvCWcnn85x4yD+PCXc7j7vLG8eM007r1gPM/mp/Gzhlvx/uR9EDr4zwXw7k1qm4H9pbXc9ub3TByUwt8WTsCo12HQ63hg4UQGp9u44+3vqW9aLO/oBnhhPrx2OfXuem5vuIF3Tv8S3RVvwJl/VK+f+Q/w+ewlzKp/DNdJi+HgV/DUHHjnRqhSr8+yWnfot4g02hb0ws3wnwvhxXPVePV5j8H/rYPJV6o1ocOR2B/O+H9w82Y1pW/9C+qd9pvHwdv+gGUIXwOseUJ15fu+gDPuhp99BYNntbtbRpKVZ3znsWfea9BQB8+dCdvfBaC4up4DZU5mDk1vtd+EgSlMy03jha8P4q0sUAeozInqxWNJ4sPvi1i2s4Tbzh4ZutgAhBD86YKxTBqUwm1vbmH/kCvUafvL/xra5r3NR/nzxzs5d0I2fzp/HPrAosmJFiPPXjmVU0ZkcsfbW9m+8St12bnM0aF9l+8u4aZXNjJ2QDIvXHNSaLaiXie45/yxnD8xh79/sjts+OVgmZObl2xiVP8k7r9wQrOY6uB0G7f/YCTLd5fy/pbCVvt6fX7+s+Ygs4amM7J/6zjmzKHpHKtxc6CsdWGq0oB77TMOfcBUUHx4j25m/cFKZoU5v6wmA5dPP4FPtxc3ZlTVlsB3/4bxP6bAPJT7P97JycMzuHhqc7Ny1axcLjxxAI98sYcVu0vUG/7FL6o3kbd+ypZDZfz8vxsYnpXIM1dOxWZWh+0Meh33LRjHwikDeWTZXjUEoiiw+mH1/E4exMGLPuKCpSb0OsGSG2Y0O7d/MmMw9y4Yy7Kdx7jtuyT8P1sNc26FLUvgienUff8eN/5nAyaDjievmNwsRdVmNnDvgnEcKnfx1Mr94DgG7/6fGpIt3wvn/oNfpT3F5+azuHBG64yuWUPTKSadLwfeCDdvUd37trfh8cl4P/4NxvpyGXLpNH6fmm+95HJ45hT1rnrWn+CXG2HKVZ1fFs2aBvMfgJ9/A4OmqaGIJ6bD1nMAoYQAABNISURBVLdUwQ57bL8aXnlyFnz6e1XAb1qrOn5Dx3fkfonqwN1B2wS4YYU6wPjmVbDsHtbsV7vDM4a0vuAArp87hPqqYlwvXKA6nyvegOQB1NQ38KcPdzBhYDLXzG6dC2w26Hly0WQSTHquX7ID1/RfqjHRA6v4fMcxfv3GFmYMSeMfP56ITtd8oMpi1PP0T6ZwyohMSnZ/S5V9SCif/cPvC7nh5fWM6G/n5WumhcQ8iE4neODiCUzPS+PXb2xplupV5fJw/cvr0ekEz/xkCglh1l28alYuEwelcO+HO1pNxvp4WzGF1fVcPSs37HcVDFmtyW/d8wo59DgU9KCINHPoAyYDULTja+oafG2eX1fNzEUnBC98fVB94pvHwOdGmXsbv39nGwrwlx+NbzWYKYTgzz8az8isRH71+mb2l9aqA7DnPgR7P2PvCzeQbjXw0jUnkRTmHLn/wvHMH9+fhz7axN5/XQzL/ghjFrDl7De4cIl6M3/1+hlhEwJ+MjOX238wknc3F3LXx/tRzvh/cP2XKLZMEt6+kl9X3ce/59nJSUlote+c4RlcMs6GbuXf8D82Gb5/A2b9EhZvZPegH/PprjKunDEYq6m1nowbkIzVpGddfoU6TnH2vbB4A0y4BP13z7HK/CvOOvJoz4V0jwNtZLlsfUuNGTsCM9Tc1eq05rm/UWPkCSldf+9+o9TQzL5l8Omd8L9r4dM/qKPng6aBrR+4a6Bwk5q9UpGvDv5cugRGzjuuzJmgiJQ43JA0GK7+UM0+Wf0PRtu/ZqDl54xpI/vijP71vG39M6aaEhqueBNjIJb94Ke7Ka918/xVJ4XcdUuykxN4/LLJXPX8t1y4bjjvWLIof/tObiz/LeNyknn2yqltTrKxGPU8vehEfH/dx3uV03n36TXoBKzNr2DyCSm8cM20Nmdcmg16nrtqKoueW8eN/93AjacMZUT/RB5ZtoeCijpevOYkBqVZw+6r1wn++qPxnPfP1fz14138baEa6mnw+Xn48z2MzErkjNFZYffNy7CRlWRmzf5yrpjefJJRSRwLusWoJ8VqpLimyQzQxP6QNJC6A98Co5mWFz57qH+yhR9OyOb17w7z82mJZH73bxi3kDcPWFi1p5R7zh/b5m+VYFJv/Bf+6xsufmoNt509kvqGU/FyITfwNvMGJ2Ozhc98Meh1PHq6mfJDf6RfyRH+nXAV39QuYvlzm8lOTuC/100Pn90V4KbThuGo9/LUyv043V7OGdefl8X9nNjwMr+yfIDxw3mwbQ6MnA+ZI9RxpMpDkL+c+w9+gtDX8a1uJlNueBx9pjph7uHPN2A3GfjpnPCTpYx6HVNz01h3oIlhSBkEC/7J1sFXk/+/uzj/8Gvw5H8g5QTIngi2TDjpesjq3VU6tSHoxd9D/gqwZahCO/R09QfqhCvuNMPOhCGnqcK+/nk1g2TdU42vCx0MmgGn/FadsKM//mnj6Ta1+3gsOOXaYFbDRNmTGPbR7bxn+A267T4Ye2FjaqDfB1vfRLf0twzQe7ms7g7mHBrAzcPVwciX1xzi6lm5jB/YfurWzKHpvHztNG55fTN/qZ3HvcYX+dWQYq5adHYrd90SS+UeUJxkjT+VmqIGhBD89pxRXDsnr8Nc20SLkZevnc5d727jn8vV/N3sZAsvXzutTbcYZExOEtednMfTK/M5dWQm88Zn86/l+zlQ5uS5K6e2eQMTQjBzSDqr95W3SpE7WlmHUS+i1iXubQamJlBQ2WIweeAUkvd8y4gsO+ntfO7FZwzn463FbFxyH2d768kf/XPuXrKdGUPS+EkHs28Hp9t488aZ3PTqJn7/jjreMi33en6cN5KUNX+F5/LhvEchZ1LjTvU1sPZJjKv/QZY5ia8mP8fSAzlUVri4ZnYevzx9eKeWCfztOSMx6gVPr8zn3c2FJJoNXHTh3RhH3Q/fPasawk9/13wnWyZi4iV8br+A6z918cvNCreeBZ9sK+aT7cXccuaIdqtxTs9L44FPd1Ph9JBma9xuv9KfWxp+wcTrx5NX9InaGz62A+rWqNd1L6MNQT/rT+q/3kanhxE/UP95PVC2G1wVajw+fThYupe7bNDryEmxcLjpzE8hKBh2Gb901/Cs/VW1h/DZnWrvQG+Cw2uh+gjkTMZw0XPkfF7Dw8v28M3+MjYcqmTSoBR+N39Up44/Y0g6K28/jT0F4/C9+QmL9W+D5dqOdzy8BoAzzl7AGV1YqSg5wchjl53Ib+eNotLpYURWYqcnXfzqjBF8d6CCxUs2Mfmbg3x7oIILJuVwxuh+7e43Y0g6724uZH+pk2H9GsdTjlS6GJCS0ObNQOsMSElgf2nzsQNfzlSydrzHmQPan1g2NNPOH07NZM7qd1hlmcsvXi8l0WLgsUtPbBWOC8eQTDsfLZ7D3pJa9DrB0EwbQsyCE8bCBzerIdKcyeq8jNpj6nnV4IIxFyDm/Z25iVnM7cJnFkLw67NHcs3sPA5XuBjezx6K1XPqHaoJcxSraYaKD5IGqEso6nScqShcWLqFx77Yy5r9ZWw5Us34AcnceOqQdo85Y4ja0/n2QHmz/P2jgZtp//4D4IQbYPoNXfhEXUcbgh4NDKY2c667wwlp1uaCDqzZX85GZQSll39GeuUKdYDl2A7wedQ4+9n3wugFoNPx94U+spIsrNhdwsVTB/G7+aOOqyaJyaBjXG4WnHwLfPJbNeMld077Ox38Ws1n7uY6ogNSEhgQJp7ZHgkmPS/+dBp/+WgnGw9X8rO5Q/j12SM7nJgSdP9r88ubCXpBhavN0EE8MDDVyqo9Zc16JntsUxgNnJ2wEzit3f2v1C0F4eEV08XMGpjOneeOCU3a6gw6nWg9UD36PMibq/Z6d32kThBKSFVTfE/8SSjO313SbKZmbjmEEJCUrf5r9ZLggYUTGZRq5ZNtxfxwQjZ3/nBMh9fU+AEpWIw61uZXNBf0qnrSbaaw40KRQAp6hDkhzcpn25vXgliTX06azcSI/imQ8yMY+6M29zcb9Px+/mh+P390m9t0iilXqZOSVv6tfUH3eWH/FzDqvMjMtA1DksXI/RdN6HjD/9/e2cdWVZ4B/PdAoRSE0tICrVCgQIfgCGhV0BF0WCF+oFkcM5tfy5xxzsxscYmLy7K4mW1ublmmyyDbMrfMjOmcsoGZH4goUqROpQMt/VBssVJaWr5aLC3P/nhP22u57T39Oqfn9vklNz338t7218Ppc9/znuc8Twyzpoxn2qRUdr1/5FPFuqobW1idG93mHImYkZFGy+l2Gk523dyytTGbLJ3EguaS3t/c0oS8sQEWrmXDupsHV2xcOlz2LfcYZoweJXy7qMB/kTLcxOii2Zln5bJXH2lmRkbfJi2DSfSyXCJOXuYEGk62cuKTNsDdBr2zsoFl+Zm+TmsHjTFpcNm9brZ0YGfP46p3ubzigquCcxsERIRl+VMormrozEc/duo0R062kpfEM/SOHPOqmGWX4vcbKR27lHEHXnEf0D2xa71LAFhx31BrJgXL506h7NDxT6WJVtSdYO7UHlKmA8ACesDMzXZ/cPsPuZsuqupPUnv0FJfNywpepvCrLlto+8M9j3lvs8sMyL88KKtB45I5Uzh8vCsfvexjt88XxMldTxY+49UXKfOOr9a2M5R80Ej9jCJornc3wsSjpQmKf+vK4eb07WxopLI8ZlkP3ITh42OnmD81vOPLAnrALMx1F1b3fnQMgNfKXc2RFfOyg5cZO8HdFFG51V187U5bK+zZCAVrItlDtOPCVUc++nu1bp8vyEnegJ6TPo6JqSmUewH9nZomWk63M3nJda5aYulT8d/46iPuTOzy+wO0jTafPTedc1JTeN1bdqmoc2Wb59sMfeRw7uQ00tPGsO8jV8vl1fLDzMxMI29KSMsAF90BE3Nh831nn46XbXGzugtuDcdtgMzJmsCMjLTOaxb7ao+RnjaG6X24yBc1RIQFORMpPeiOr21ldYweJVw8PxcWXQ97n4aT3QqXNVS6FN0lX/50WqHRKymjR7EsP5Pt+12Zibc/bALgvNzwKnlaQA8YEWHJzMnsrGzgaMtptu+vp+i86eEJpZ4DVz8Mh0phx6+6Xj/T7mp5ZMyBuavC8xsAIsI1i3PYUVFP48lWXquo56LZGX0u7xs1luVPYU/NUY62nOa5/33MsvxMl1N96b2u5MTrv+ka3N7mbntPSXP1zY0+cdXC6dQ0tlB68Cg7qxrIyxzf50yuwWRAAV1E1ohImYhUiIidq/mkaOE0Pmho5sF/7aO1/Qw3LPXZYm6oWHCtu1lq60NdLcu2/tgF+VU/8F9OYRhy3eJc2s4oD215l+ojLVyxoPf89WTgc/OyaD+jPLR5H1WHT3a1MMwugMXrXBerD3a4JbVnvwnVxXDNL2BSyMdhBFm9aDqpKaNY/0oVxZUNcWvlBEm//1JFZDTwGFAE1AC7RWSTqgZfwCBirDnfle/8x39rWFmQzeIZAyhdMBiIwNpHXUGmZ+5yfUmbG1yOcC8plFFgUe4kVszP4qk3axg/djSrF4V4NhQQhbMzmZM1gb+X1JCTPo4blnbV12fNT10D8cevhdRJcKrJlZdevC484QiTPn4Mty6f1dkl64uFM0P1kUT1qnt8o8hy4Iequtp7/j0AVf1JT+8pLCzUkpIEubAjhNcr6nnx3TruWpnfpxs3hpT2Nlf24OCbrvDY0lt67U4UFeqOnWL99ipWLZjKpWFkE4VA+aHj/HXXh6wrnNl5Ib6T5iMuo+V4rfvAnndlOJJJQktrO4+9XEF+9gS+cMHQtGcUkTdVtTDhuAEE9BuBNap6h/f8FuASVb2n27g7gTsB8vLyLjxw4EC/fp5hGMZIxW9AH8j0K96VpbM+HVR1g6oWqmphdnYIqXmGYRgjhIEE9BogdsFoBnB2NwLDMAwjEAYS0HcD80VkjoiMBW4CNg2OlmEYhtFX+p3loqptInIP8B9gNPBHVd2b4G2GYRjGEDGgBGNV3QJsGSQXwzAMYwBEPyfNMAzDACygG4ZhJA0W0A3DMJKEft9Y1K8fJnIYGG53FmUB9QlHDQ+i5ArR8o2SK0TLN0quMDx9Z6lqwht5Ag3owxERKfFzB9ZwIEquEC3fKLlCtHyj5ArR843FllwMwzCSBAvohmEYSYIFdNgQtkAfiJIrRMs3Sq4QLd8ouUL0fDsZ8WvohmEYyYLN0A3DMJKEERfQRSRTRF4QkXLva0YP4/JE5HkReVdE9onI7GBN/bt6YyeJyEEReTRIx24OCX1FZImI7BSRvSKyR0S+FLBjr20TRSRVRDZ6/74rjP/3bj6JfL/jHZ97ROQlEZkVhqfn4qslpYjcKCIqIqFlkvhxFZF13r7dKyJPBO3YL1R1RD2Ah4H7ve37gZ/1MG4bUORtnwOMH66u3r//GngCeHQ471ugAJjvbecCtcDkgPxGA5VAPjAWeAdY2G3M3cDvvO2bgI0h7k8/vld0HJvAN8Ly9ePqjZsIbAeKgcLh6grMB94CMrznU8M6DvryGHEzdOB64HFv+3Hghu4DRGQhkKKqLwCo6glVbQ5OsZOErgAiciEwDXg+IK+eSOirqvtVtdzb/gioA4LqfHIxUKGqVaraCvwN5xxL7O/wFLBKROI1cwmChL6q+nLMsVmM60sQBn72LcCPcB/8p4KU64Yf168Dj6lqI4Cq1gXs2C9GYkCfpqq1AN7XeG3gC4AmEXlaRN4SkZ97TbGDJqGriIwCHgG+G7BbPPzs205E5GLcDKkyADeAc4HqmOc13mtxx6hqG3AUCKuVux/fWL4GPDekRj2T0FVElgIzVfXfQYrFwc9+LQAKRGSHiBSLyJrA7AbAgMrnDldE5EUgXnv3B3x+ixRgBbAU+BDYCNwO/GEw/GIZBNe7gS2qWh3ERHIQfDu+Tw7wF+A2VT0zGG5+fmyc17qneflqrRgQvl1E5GagEFg5pEY906urN/H4Fe7vKGz87NcU3LLL5biznldF5HxVbRpitwGRlAFdVXtsYy4ih0QkR1VrvaAS71SqBnhLVau89zwDLGMIAvoguC4HVojI3bi1/rEickJVe7woFbIvIjIJ2Ax8X1WLh8KzB/y0TewYUyMiKUA6cCQYvbPw1eZRRK7EfaCuVNVPAnLrTiLXicD5wDZv4jEd2CQia1W1JDBLh9/joFhVTwPvi0gZLsDvDkaxf4zEJZdNwG3e9m3As3HG7AYyRKRjbffzwL4A3LqT0FVVv6Kqeao6G7gP+PNQBXMfJPT12hX+E+f5ZIBu4K9tYuzvcCOwVb2rYiGQ0NdbxlgPrA15nbdXV1U9qqpZqjrbO1aLcc5BB/OErh7P4C44IyJZuCWYqkAt+0PYV2WDfuDWQ18Cyr2vmd7rhcDvY8YVAXuAUuBPwNjh6hoz/nbCzXJJ6AvcDJwG3o55LAnQ8WpgP27d/gHvtQdxwQVgHPAkUAG8AeSHfLwm8n0ROBSzLzcNV9duY7cRUpaLz/0qwC9xE7lS4KYwjwO/D7tT1DAMI0kYiUsuhmEYSYkFdMMwjCTBArphGEaSYAHdMAwjSbCAbhiGkSRYQDcMw0gSLKAbhmEkCRbQDcMwkoT/A9KueXVymnFRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "np.set_printoptions(precision = 4)\n",
    "def grating_fourier_harmonics(order, fill_factor, n_ridge, n_groove):\n",
    "    \"\"\" function comes from analytic solution of a step function in a finite unit cell\"\"\"\n",
    "    #n_ridge = index of refraction of ridge (should be dielectric)\n",
    "    #n_ridge = index of refraction of groove (air)\n",
    "    #n_ridge has fill_factor\n",
    "    #n_groove has (1-fill_factor)\n",
    "    # there is no lattice constant here, so it implicitly assumes that the lattice constant is 1...which is not good\n",
    "\n",
    "    if(order == 0):\n",
    "        return n_ridge**2*fill_factor + n_groove**2*(1-fill_factor);\n",
    "    else:\n",
    "        #should it be 1-fill_factor or fill_factor?, should be fill_factor\n",
    "        return(n_ridge**2 - n_groove**2)*np.sin(np.pi*order*(fill_factor))/(np.pi*order);\n",
    "\n",
    "def grating_fourier_array(num_ord, fill_factor, n_ridge, n_groove):\n",
    "    \"\"\" what is a convolution in 1D \"\"\"\n",
    "    fourier_comps = list();\n",
    "    for i in range(-num_ord, num_ord+1):\n",
    "        fourier_comps.append(grating_fourier_harmonics(i, fill_factor, n_ridge, n_groove));\n",
    "    return fourier_comps;\n",
    "\n",
    "def fourier_reconstruction(x, period, num_ord, n_ridge, n_groove, fill_factor = 0.5):\n",
    "    index = np.arange(-num_ord, num_ord+1);\n",
    "    f = 0;\n",
    "    for n in index:\n",
    "        coef = grating_fourier_harmonics(n, fill_factor, n_ridge, n_groove);\n",
    "        f+= coef*np.exp(cmath.sqrt(-1)*np.pi*n*x/period);\n",
    "        #f+=coef*np.cos(np.pi*n*x/period)\n",
    "    return f;\n",
    "\n",
    "def fourier_reconstruction_general(x, period, num_ord, coefs):\n",
    "    '''\n",
    "    overloading odesn't work in python...fun fact, since it is dynamically typed (vs statically typed)\n",
    "    :param x:\n",
    "    :param period:\n",
    "    :param num_ord:\n",
    "    :param coefs:\n",
    "    :return:\n",
    "    '''\n",
    "    index = np.arange(-num_ord, num_ord+1);\n",
    "    f = 0; center = int(len(coefs)/2); #no offset\n",
    "    for n in index:\n",
    "        coef = coefs[center+n];\n",
    "        f+= coef*np.exp(cmath.sqrt(-1)*2*np.pi*n*x/period);\n",
    "    return f;\n",
    "\n",
    "def grating_fft(eps_r):\n",
    "    assert len(eps_r.shape) == 2\n",
    "    assert eps_r.shape[1] == 1;\n",
    "    #eps_r: discrete 1D grid of the epsilon profile of the structure\n",
    "    fourier_comp = np.fft.fftshift(np.fft.fft(eps_r, axis = 0)/eps_r.shape[0]);\n",
    "    #ortho norm in fft will do a 1/sqrt(n) scaling\n",
    "    return np.squeeze(fourier_comp);\n",
    "\n",
    "# plt.plot(x, np.real(fourier_reconstruction(x, period, 1000, 1,np.sqrt(12), fill_factor = 0.1)));\n",
    "# plt.title('check that the analytic fourier series works')\n",
    "# #'note that the lattice constant tells you the length of the ridge'\n",
    "# plt.show()\n",
    "\n",
    "L0 = 1e-6;\n",
    "e0 = 8.854e-12;\n",
    "mu0 = 4*np.pi*1e-8;\n",
    "fill_factor = 0.3; # 50% of the unit cell is the ridge material\n",
    "\n",
    "\n",
    "num_ord = 10; #INCREASING NUMBER OF ORDERS SEEMS TO CAUSE THIS THING TO FAIL, to many orders induce evanescence...particularly\n",
    "               # when there is a small fill factor\n",
    "PQ = 2*num_ord+1;\n",
    "indices = np.arange(-num_ord, num_ord+1)\n",
    "\n",
    "n_ridge = 3.48; #3.48;              # ridge\n",
    "n_groove = 1;                # groove (unit-less)\n",
    "lattice_constant = 0.7;  # SI units\n",
    "# we need to be careful about what lattice constant means\n",
    "# in the gaylord paper, lattice constant exactly means (0, L) is one unit cell\n",
    "\n",
    "\n",
    "d = 0.46;               # thickness, SI units\n",
    "\n",
    "Nx = 2*256;\n",
    "eps_r = n_groove**2*np.ones((2*Nx, 1)); #put in a lot of points in eps_r\n",
    "border = int(2*Nx*fill_factor);\n",
    "eps_r[0:border] = n_ridge**2;\n",
    "fft_fourier_array = grating_fft(eps_r);\n",
    "x = np.linspace(-lattice_constant,lattice_constant,1000);\n",
    "period = lattice_constant;\n",
    "fft_reconstruct = fourier_reconstruction_general(x, period, num_ord, fft_fourier_array);\n",
    "\n",
    "fourier_array_analytic = grating_fourier_array(Nx, fill_factor, n_ridge, n_groove);\n",
    "analytic_reconstruct = fourier_reconstruction(x, period, num_ord, n_ridge, n_groove, fill_factor)\n",
    "\n",
    "\n",
    "plt.figure();\n",
    "plt.plot(np.real(fft_fourier_array[Nx-20:Nx+20]), linewidth=2)\n",
    "plt.plot(np.real(fourier_array_analytic[Nx-20:Nx+20]));\n",
    "plt.legend(('fft', 'analytic'))\n",
    "plt.show()\n",
    "\n",
    "plt.figure();\n",
    "plt.plot(x,fft_reconstruct)\n",
    "plt.plot(x,analytic_reconstruct);\n",
    "plt.legend(['fft', 'analytic'])\n",
    "plt.show()\n",
    "theta_inc = 0;\n",
    "\n",
    "\n",
    "wavelength_scan = np.linspace(0.5,2.3,300)\n",
    "## construct permittivity harmonic components E\n",
    "#fill factor = 0 is complete dielectric, 1 is air\n",
    "\n",
    "\n",
    "##construct convolution matrix\n",
    "E_conv = np.zeros((2 * num_ord + 1, 2 * num_ord + 1));\n",
    "E_conv = E_conv.astype('complex')\n",
    "p0 = Nx;\n",
    "p_index = np.arange(-num_ord, num_ord + 1);\n",
    "fourier_array = fft_fourier_array; #_analytic;\n",
    "for prow in range(2 * num_ord + 1):\n",
    "    # first term locates z plane, 2nd locates y coumn, prow locates x\n",
    "    for pcol in range(2 * num_ord + 1):\n",
    "        pfft = p_index[prow] - p_index[pcol];\n",
    "        E_conv[prow, pcol] = fourier_array[p0 + pfft];  # fill conv matrix from top left to top right\n",
    "\n",
    "## FFT of 1/e;\n",
    "inv_fft_fourier_array = grating_fft(1/eps_r);\n",
    "##construct convolution matrix\n",
    "E_conv_inv = np.zeros((2 * num_ord + 1, 2 * num_ord + 1));\n",
    "E_conv_inv = E_conv_inv.astype('complex')\n",
    "p0 = Nx;\n",
    "p_index = np.arange(-num_ord, num_ord + 1);\n",
    "for prow in range(2 * num_ord + 1):\n",
    "    # first term locates z plane, 2nd locates y coumn, prow locates x\n",
    "    for pcol in range(2 * num_ord + 1):\n",
    "        pfft = p_index[prow] - p_index[pcol];\n",
    "        E_conv_inv[prow, pcol] = inv_fft_fourier_array[p0 + pfft];  # fill conv matrix from top left to top right\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'bEr' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-3-6084a00e4fd4>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     27\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     28\u001b[0m     \u001b[1;31m## ======================== off diagonal terms =================================================== ##\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 29\u001b[1;33m     \u001b[0mbE\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0minv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbEr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m;\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     30\u001b[0m     \u001b[0mG\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mj\u001b[0m\u001b[1;33m*\u001b[0m \u001b[0mEzxzz\u001b[0m \u001b[1;33m@\u001b[0m \u001b[0mKX\u001b[0m\u001b[1;31m#j*bslash(Ezz,Ezx) @ KX;\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     31\u001b[0m     \u001b[1;31m#G = j*(Ezx/Ezz)@KX #we should not do pointwise division of these epsilon matrices.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'bEr' is not defined"
     ]
    }
   ],
   "source": [
    "I = np.identity(2*num_ord+1);\n",
    "spectra = list();\n",
    "spectra_T = list();\n",
    "# E is now the convolution of fourier amplitudes\n",
    "for lam0 in wavelength_scan:\n",
    "    k0 = 2*np.pi/lam0; #free space wavelength in SI units\n",
    "    #print('wavelength: ' + str(lam0));\n",
    "    ## =====================STRUCTURE======================##\n",
    "\n",
    "    ## Region I\n",
    "    n1 = 1;#cmath.sqrt(-1)*1e-12; #apparently small complex perturbations are bad in Region 1, these shouldn't be necessary\n",
    "    ## Region 2; transmission\n",
    "    n2 = 1;\n",
    "\n",
    "    #from the kx_components given the indices and wvln\n",
    "    #2 * np.pi * np.arange(-N_p, N_p + 1) / (k0 * a_x)\n",
    "    indices = np.arange(-num_ord, num_ord + 1)\n",
    "    kx_array = (n1*np.sin(theta_inc) + indices*(lam0 / lattice_constant)); #0 is one of them, k0*lam0 = 2*pi\n",
    "\n",
    "    ## IMPLEMENT SCALING: these are the fourier orders of the x-direction decomposition.\n",
    "    KX = np.diag(kx_array); #singular since we have a n=0, m= 0 order and incidence is normal\n",
    "\n",
    "    ## construct matrix of Gamma^2 ('constant' term in ODE):\n",
    "    #use fast fourier factorization rules here...\n",
    "    Q = I- KX @ bslash(E_conv, KX);\n",
    "    PQ = -bslash(E_conv_inv, Q);\n",
    "    \n",
    "    ## ======================== off diagonal terms =================================================== ##\n",
    "    bE = np.linalg.inv(bEr);\n",
    "    G = j* Ezxzz @ KX#j*bslash(Ezz,Ezx) @ KX;\n",
    "    #G = j*(Ezx/Ezz)@KX #we should not do pointwise division of these epsilon matrices.\n",
    "    H = -j*KX @Exzzz; #j*KX @bslash(Ezz, Exz);\n",
    "    M = np.linalg.inv(bE);\n",
    "    K = -(B + H@np.linalg.inv(bE)@G);\n",
    "    C = -np.linalg.inv(bE)@G - H@np.linalg.inv(bE);\n",
    "    Z = np.zeros_like(M);\n",
    "    I = np.eye(M.shape[0], M.shape[1]);\n",
    "    \n",
    "    ## ================================================================================================##\n",
    "    eigenvals, W = LA.eig(PQ); #A is NOT HERMITIAN\n",
    "    #we should be gauranteed that all eigenvals are REAL\n",
    "    eigenvals = eigenvals.astype('complex');\n",
    "    lambda_matrix = np.diag((np.sqrt((eigenvals))));\n",
    "\n",
    "    ## THIS NEGATIVE SIGN IS CRUCIAL, but I'm not sure why\n",
    "    V = -em.eigen_V(Q, W, lambda_matrix)\n",
    "    kz_inc = n1;\n",
    "    ## ================================================================================================##\n",
    "\n",
    "    ## scattering matrix needed for 'gap medium'\n",
    "    #if calculations shift with changing selection of gap media, this is BAD; it should not shift with choice of gap\n",
    "    Wg,Vg, Kzg = hl.homogeneous_1D(KX, 1, m_r = 1)\n",
    "    ## reflection medium\n",
    "    Wr,Vr, Kzr = hl.homogeneous_1D(KX, 1, m_r = 1)\n",
    "    ## transmission medium;\n",
    "    Wt,Vt, Kzt = hl.homogeneous_1D(KX, 1, m_r = 1)\n",
    "\n",
    "    ## S matrices for the reflection region\n",
    "    #Ar, Br = sm.A_B_matrices(Wg, Wr, Vg, Vr);\n",
    "    Ar, Br = sm.A_B_matrices_half_space(Wr, Wg, Vr, Vg);  # make sure this order is right\n",
    "\n",
    "    S_ref, Sr_dict = sm.S_R(Ar, Br);  # scatter matrix for the reflection region    ## calculating A and B matrices for scattering matrix\n",
    "    Sg = Sr_dict;\n",
    "\n",
    "    ## define S matrix for the GRATING REGION\n",
    "    A, B = sm.A_B_matrices(W, Wg,  V, Vg);\n",
    "    S, S_dict = sm.S_layer(A, B, d, k0, lambda_matrix)\n",
    "    Sg_matrix, Sg = rs.RedhefferStar(Sg, S_dict)\n",
    "\n",
    "    ## define S matrices for the Transmission region\n",
    "    At, Bt = sm.A_B_matrices_half_space(Wt, Wg, Vt, Vg);  # make sure this order is right\n",
    "    St, St_dict = sm.S_T(At, Bt); #scatter matrix for the reflection region\n",
    "    Sg_matrix, Sg = rs.RedhefferStar(Sg, St_dict)\n",
    "\n",
    "    #check scattering matrix is unitary\n",
    "    #print(np.linalg.norm(np.linalg.inv(Sg_matrix)@Sg_matrix - np.matrix(np.eye(2*(2*num_ord+1)))))\n",
    "\n",
    "    ## ======================== CALCULATE R AND T ===============================##\n",
    "    K_inc_vector =  n1*np.matrix([np.sin(theta_inc), \\\n",
    "                                         0, np.cos(theta_inc)]);\n",
    "    #K_inc isn't even used for anyting...\n",
    "\n",
    "    #cinc is the incidence vector\n",
    "    cinc = np.zeros((2*num_ord+1, )); #only need one set...\n",
    "    cinc[num_ord] = 1;\n",
    "    cinc = cinc.T;\n",
    "    cinc = np.linalg.inv(Wr) @ cinc;\n",
    "    ## COMPUTE FIELDS: similar idea but more complex for RCWA since you have individual modes each contributing\n",
    "    reflected = Wr @ Sg['S11'] @ cinc;\n",
    "    transmitted = Wt @ Sg['S21'] @ cinc;\n",
    "\n",
    "    ## reflected is already ry or Ey\n",
    "    rsq = np.square(np.abs(reflected));\n",
    "    tsq = np.square(np.abs(transmitted));\n",
    "\n",
    "    ## compute final reflectivity\n",
    "    Rdiff = np.real(Kzr)@rsq/np.real(kz_inc); #real because we only want propagating components\n",
    "    Tdiff = np.real(Kzt)@tsq/np.real(kz_inc)\n",
    "    R = np.sum(Rdiff);\n",
    "    T = np.sum(Tdiff);\n",
    "\n",
    "    #print(R);\n",
    "    spectra.append(R); #spectra_T.append(T);\n",
    "    spectra_T.append(T)\n",
    "\n",
    "plt.figure();\n",
    "plt.plot(wavelength_scan, spectra);\n",
    "plt.plot(wavelength_scan, spectra_T)\n",
    "plt.plot(wavelength_scan, np.array(spectra)+np.array(spectra_T))\n",
    "plt.legend(['reflection', 'transmission'])\n",
    "plt.show()\n",
    "\n",
    "\n",
    "\n",
    "# ## now we can form the linear system to solve for the amplitude coeffs\n",
    "# WV1 = np.block([[W, np.matmul(W,X)],[V, -np.matmul(V,X)]]);\n",
    "#\n",
    "# WV2 = np.block([[np.matmul(W,X), W],[np.matmul(V,X), -V]])\n",
    "#\n",
    "# print(WV1.shape)\n",
    "# '''\n",
    "# the conditioning of WV1 and WV2 cannot suck because we're going to invert one of them.\n",
    "# '''\n",
    "# print('condition of WV1: '+str(np.linalg.cond(WV1))) #condition number of this is HUGE!!! caused by X\n",
    "# print('condition of WV2: '+str(np.linalg.cond(WV2))) #condition number of this is HUGE!!!\n",
    "#\n",
    "# F1 = np.linalg.inv(WV1);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
